Analyses for stickleback landscape parasitology study
Pascal I. Hablützel and Io S. Deflem
2022-11-27
1. Description
This R Markdown document describes the analyses performed for the manuscript entitled “Environmental pollution correlates with parasite infection across a riverine landscape” by Io S. Deflem, Seppe Marchand, Federico C.F. Calboli, Joost A.M. Raeymaekers, Filip A.M. Volckaert and Pascal I. Hablützel.
The analyses were run in R 4.2.2
2. Study area and sampling
Up to thirty 0+ three-spined sticklebacks were sampled at 37 locations in the Dijle and Demer basins in Belgium during autumn 2016 under a permit of the Flemish Agency Nature and Forest. Both basins together cover a continuous surface area of 3,624 km² with the furthest two sampling sites being located 117 km apart (distance measured along rivers). All locations included small and relatively slow flowing streams (drop off from highest to lowest point is 18 m) covering a wide range of ecological, hydromorphological, and physico-chemical characteristics. Fish were caught using a dip net.
3. Setting up working environment
# Empty environment
rm(list=ls())
# Set working directory to location where script is stored
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) # requires installation of package "rstudioapi"
# Loading required libraries
require(BAS)
require(boral)
require(car)
require(corrplot)
require(ggplot2)
require(gplots)
require(vegan)4 Loading and preparing host and parasite data
Fish were euthanized with a lethal dose of MS222 on the day of capture, following directions of the KU Leuven Animal Ethics Commission, and stored at -20 °C. Fish were kept in separate containers per site at all times. Lab based parasite screening of thawed fish involved placing individual fish in 5 or 10 ml cryo-tubes with 1 or 2 ml of distilled water. Following a vigorous shake of 10 s, the liquid was poured into a Petri dish and ectoparasites were identified and counted using a stereomicroscope. Fish were rinsed and checked again for the presence of ectoparasites on skin and fins. The intestines were examined for endoparasites. Before dissection, fish weight (± 0.1 mg) and standard length (± 1 mm) were recorded. To quantify body condition, we calculated the scaled mass index (SMI; Maceda-Veiga et al., 2014; Peig & Green, 2009). Sex was determined during dissection by inspection of gonad development. A total of 668 fish were dissected, which amounts to approximately 20 fish per location, with the exception of seven locations where only 10 fish were screened for the presence of macroparasites. Ecto- and endoparasites were morphologically identified to species level whenever possible.
# Parasite data
data <- read.csv("data_2016_2303.csv", sep=';')
data$site <- as.factor(data$site)
# Calculate parasite parameters
#names(data)
#parasite data is overdispersed (mostly so for Trichodina), if using average abundance data, species matrix needs to be transformed
datao <- na.omit(data[,c(1,22:24,26:32)]) #remove fish without parasite counts
ddata <- dispweight(datao[,-1]) #correct for overdispersion of the parasite count data
avab <- aggregate(ddata, by = list(datao[,1]), function(x){mean(x, na.rm =T)})
prev = aggregate(data[,c(22:24,26:32)], by = list(data[,1]), function(x){sum(x >0, na.rm = T)/length(x)})
medin = aggregate(data[,c(22:24,26:32)], by = list(data[,1]), function(x){median(x[x >0], na.rm = T)})
pa = aggregate(data[,c(22:24,26:32)], by = list(data[,1]), function(x){ifelse(mean(x, na.rm =T)>0,1,0)})
avab[is.na(avab)] <- 0
prev[is.na(prev)] <- 0
medin[is.na(medin)] <- 0
# Host condition data
avcondition <- aggregate(data$SMI, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2]
avlength <- aggregate(data$length, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2]
# Parasite index
sgyr <- 1:nrow(data)
stri <- 1:nrow(data)
sglu <- 1:nrow(data)
scon <- 1:nrow(data)
scysl <- 1:nrow(data)
spro <- 1:nrow(data)
saca <- 1:nrow(data)
scam <- 1:nrow(data)
sang <- 1:nrow(data)
for(j in 1:nrow(data)){
sgyr[j] <- data$Gyr[j]/sd(data$Gyr, na.rm=T)
stri[j] <- data$Tri[j]/sd(data$Tri, na.rm=T)
sglu[j] <- data$Glu[j]/sd(data$Glu, na.rm=T)
scon[j] <- data$Con[j]/sd(data$Con, na.rm=T)
scysl[j] <- data$CysL[j]/sd(data$CysL, na.rm=T)
spro[j] <- data$Pro[j]/sd(data$Pro, na.rm=T)
saca[j] <- data$Aca[j]/sd(data$Aca, na.rm=T)
scam[j] <- data$Cam[j]/sd(data$Cam, na.rm=T)
sang[j] <- data$Ang[j]/sd(data$Ang, na.rm=T)
}
PI <- 1:nrow(data)
for(j in 1:nrow(data)){
PI[j] <- 10/max(sgyr, na.rm=T)*sgyr[j] + 10/max(stri, na.rm=T)*stri[j] + 10/max(sglu, na.rm=T)*sglu[j] +
10/max(scon, na.rm=T)*scon[j] + 10/max(scysl, na.rm=T)*scysl[j] + 10/max(spro, na.rm=T)*spro[j] +
10/max(saca, na.rm=T)*saca[j] + 10/max(scam, na.rm=T)*scam[j] + 10/max(sang, na.rm=T)*sang[j]
}
PI_ecto <- 1:nrow(data)
for(j in 1:nrow(data)){
PI_ecto[j] <- 10/max(sgyr, na.rm=T)*sgyr[j] + 10/max(stri, na.rm=T)*stri[j] + 10/max(sglu, na.rm=T)*sglu[j]
}
PI_endo <- 1:nrow(data)
for(j in 1:nrow(data)){
PI_endo[j] <- 10/max(scon, na.rm=T)*scon[j] + 10/max(scysl, na.rm=T)*scysl[j] + 10/max(spro, na.rm=T)*spro[j] +
10/max(saca, na.rm=T)*saca[j] + 10/max(scam, na.rm=T)*scam[j] + 10/max(sang, na.rm=T)*sang[j]
}
avPI <- aggregate(PI, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2]
avPI_ecto <- aggregate(PI_ecto, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2]
avPI_endo <- aggregate(PI_endo, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2]5 Loading and preparing environmental and spatial data
Physico-chemical data was provided by the Flemish Environmental Agency (VMM). Each fish sampling site was chosen at or near an environmental monitoring site of VMM. Water parameters include water temperature, pH, conductivity, dissolved oxygen (O2), saturation with dissolved oxygen, and Biochemical and Chemical Oxygen Demand (BOD and COD). Nutrient related water parameters include measurements of nitrate (NO3-), nitrite (NO2-), Kjeldahl nitrogen (KjN), total nitrogen (Nt), ammonium (NH4+), and total phosphorus (Pt). Following removal of strong collinear variables (significant correlation with P < 0.05 and Pearson correlation coefficient > |0.6|; Dormann et al., 2013), six environmental physico-chemical variables were retained (temperature, conductivity, COD, saturation with dissolved oxygen, ammonium, and total nitrogen), representing different aspects of water quality. For each water parameter, the average value of the year before sampling was calculated based on monthly monitoring data. Additionally, two hydromorphological variables were included: Tthe presence of a pool-riffle pattern and meanders were noted during field sampling and these parameters were included as binary variables (presence/absence) for a representation of abiotic habitat structure. Spatial (waterway) distances were calculated using the Network Analyst toolbox in ArcGIS. Upstream distance was defined as the maximal upstream distance from the sampling location. Network peripherality was calculated as the average waterway distance of a sampling location to all other locations. Hence, a total of eight environmental and two spatial variables were included in the statistical analysis.
# Environmental data (VMM)
environment <- read.csv("Environment_update.csv", sep=';')
spavar <- read.csv("space2.csv", sep=';') #spatial variables: network centrality and upstream distance
plot(spavar$netcen); plot(density(spavar$netcen))
plot(spavar$updist); plot(density(spavar$updist))
# Environmental data (from field observations)
field_data <- read.csv("field_data.csv", sep=',')
environment2 <- cbind(environment[,c(1,49,52:53,55,57,60,63)], field_data[-c(8,10,25,27,37),33:34], spavar[,c(2,3)])
environment2$pool_riffle <- as.factor(environment2$pool_riffle)
environment2$meander <- as.factor(environment2$meander)
netcen <- spavar$netcen
updist <- spavar$updistWe used univariate generalized linear models to investigate how landscape-level effects modify infection patterns of individual parasite taxa, host size and condition. We kept the statistical models linear (as opposed to polynomial) and only considered main effects (i.e. no interaction terms) because we had no prior information from this study system that more complex models were necessary and because the study design with (only) 37 sampling sites was not intended for non-linear models. Ten explanatory variables (temperature, conductivity, COD, saturation with dissolved oxygen, ammonium, total nitrogen, the presence of pool-riffle patterns and meanders, upstream distance, and network peripherality) were included.
6. Univeriate analysis using Bayesian Model Averaging
Univariate analyses - We used generalized linear models in a BMA approach to understand how infection with individual parasite taxa relate to host characteristics (length and condition), environmental conditions as well as spatial distance among sampling sites. Parasite infection was calculated in three ways at the host population level: average abundance (mean parasites per host), prevalence (percentage of infected hosts) and median infection intensity (median number of parasites in infected hosts). We calculated the individual parasitation index (IPI) following Kalbe et al. (2002) as a measurement for total parasite abundance and species richness for each individual fish. This index was calculated for all parasite species combined, and for ecto- and endoparasite species separately. For these models, we assumed a normal error distribution (which appeared to be justified, see Supplementary Figures S1-S2) and applied a Jeffrey-Zellner-Siow prior. Model assumptions (homoscedasticity of the variances and normal distribution of the errors) were assessed using the generic model plot function in R and did onlynot clearly deviate in any of the models for rare parasites. We followed a normal distribution, and not Poisson or negative binomial, for the parasite data for the common species (Trichodina sp. and Gyrodactylus spp.) and the individual parasitation index as the parameters used are deviates from count data. Rare parasites (Glugea, Contracaecum, Anguillicoloides, and unidentified cysts) were excluded from the univariate analysis because there was not enough data to obtain a good fit of the models.For rare parasites (Contracaecum sp. and Anguillicoloides crassus), we used population-level presence-absence data assuming a binomial error distribution and a uniformly distributed BIC prior. Due to low prevalences, the other parasites were not included in the species-specific models. Explanatory variables were considered important when they had a posterior inclusion probability (PIP) of 0.5. To account for overdispersion in the parasite counts, we transformed the data by downweighting overdispersed taxa following Clarke et al. (2006) using the dispweight function in the R package vegan v2.5.6 (Oksanen et al., 2013).
6.0.1 Figure 2
# Make a matrix for PIP (Posterior Inclusion Probability)
PIP <- matrix(nrow=12, ncol=14)
rownames(PIP) <- c("Host condition", "Host length", "Temperature", "Oxygen saturation", "Conductivity", "COD", "Ammonium", "Total nitrogen", "Pool riffle pattern", "Meander", "Network peripherality", "Upstream distance")
colnames(PIP) <- c("Host condition", "Host size", "Gyrodactylus abundance", "Gyrodactylus prevalence", "Gyrodactylus infection intensity", "Trichodina abundance", "Trichodina prevalence", "Trichodina infection intensity", "Glugea", "Contracaecum", "Aguillicola",
"PI", "PI ecto", "PI endo") #Condition
bas.model <- bas.lm(avcondition ~ T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(3:12),1] <- pip[2:11,1]*sign(coef.model$postmean[2:11])
#Length
bas.model <- bas.lm(avlength ~ T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(3:12),2] <- pip[2:11,1]*sign(coef.model$postmean[2:11])
#Gyrodactylus abundance
bas.model <- bas.lm(avab$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#Gyrodactylus prevalence
bas.model <- bas.lm(prev$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),4] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#Gyrodactylus infection intensity
bas.model <- bas.lm(medin$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),5] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#Trichodina abundance
bas.model <- bas.lm(avab$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),6] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#Trichodina prevalence
bas.model <- bas.lm(prev$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),7] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#Trichodina infection intensity
bas.model <- bas.lm(medin$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),8] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#Glugea
bas.model <- bas.glm(pa$Glu ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, betaprior=g.prior(100), family=binomial)
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),9] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#Contracaecum
bas.model <- bas.glm(pa$Con ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, betaprior=g.prior(100), family=binomial)
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),10] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#Anguillicola
bas.model <- bas.glm(pa$Ang ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, betaprior=g.prior(100), family=binomial)
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),11] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#PI
bas.model <- bas.lm(avPI ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),12] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#PI ecto
bas.model <- bas.lm(avPI_ecto ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),13] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
#PI endo
bas.model <- bas.lm(avPI_endo ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),14] <- pip[2:13,1]*sign(coef.model$postmean[2:13])x = round(PIP, digits=2)
x[abs(PIP)<0.5] <- ""
x[abs(PIP)>0.5] <- "+"
heatmap.2(PIP[,-c(9,10,11)],
cellnote = x[,-c(9,10,11)],
#main = "Correlation",
notecex=1,
notecol="white",
density.info="none",
trace="none",
margins =c(10,8),
col=redblue(256),
dendrogram="both",
cexRow = 0.7,
cexCol = 0.7,
key.title = "PIP",
lhei = c(1,3),
lwid = c(0.5, 0.5),
srtCol = 45,
labCol = c("Host condition", "Host length", expression(paste(italic("Gyrodactylus"), " spp. - abundance")), expression(paste(italic("Gyrodactylus"), " spp. - prevalence")), expression(paste(italic("Gyrodactylus"), " spp. - infection intensity")), expression(paste(italic("Trichodina"), " sp. - abundance")), expression(paste(italic("Trichodina"), " sp. - prevalence")), expression(paste(italic("Trichodina"), " sp. - infection intensity")), expression("I"[PI]*" all parasites"), expression("I"[PI]*" ectoparasites"), expression("I"[PI]*" endoparasites")) ,
#Colv="NA"
) 
pdf("Figure2.pdf", width = 7.29, height = 4.5)
heatmap.2(PIP[,-c(9,10,11)],
cellnote = x[,-c(9,10,11)],
#main = "Correlation",
notecex=1,
notecol="white",
density.info="none",
trace="none",
margins =c(10,8),
col=redblue(256),
dendrogram="both",
cexRow = 0.7,
cexCol = 0.7,
key.title = "PIP",
lhei = c(1,3),
lwid = c(0.5, 0.5),
srtCol = 45,
labCol = c("Host condition", "Host length", expression(paste(italic("Gyrodactylus"), " spp. - abundance")), expression(paste(italic("Gyrodactylus"), " spp. - prevalence")), expression(paste(italic("Gyrodactylus"), " spp. - infection intensity")), expression(paste(italic("Trichodina"), " sp. - abundance")), expression(paste(italic("Trichodina"), " sp. - prevalence")), expression(paste(italic("Trichodina"), " sp. - infection intensity")), expression("I"[PI]*" all parasites"), expression("I"[PI]*" ectoparasites"), expression("I"[PI]*" endoparasites")) ,
#Colv="NA"
)
dev.off()## png
## 26.1 Variation in host condition
bas.model <- bas.lm(avcondition ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av +
pool_riffle + meander + netcen + updist,
data=environment2, prior="JZS")
plot(bas.model)
summary(bas.model)## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.00000000 1.0000 1.0000000 1.0000000 1.0000000 1.0000000
## T_av 0.12386945 0.0000 1.0000000 0.0000000 0.0000000 0.0000000
## O2_sat_av 0.05073157 0.0000 0.0000000 0.0000000 0.0000000 0.0000000
## Con_av 0.03919858 0.0000 0.0000000 0.0000000 0.0000000 0.0000000
## COD_av 0.05197281 0.0000 0.0000000 0.0000000 0.0000000 1.0000000
## NH4._av 0.04846220 0.0000 0.0000000 0.0000000 0.0000000 0.0000000
## Nt_av 0.05248386 0.0000 0.0000000 0.0000000 0.0000000 0.0000000
## pool_riffle1 0.06729856 0.0000 0.0000000 0.0000000 1.0000000 0.0000000
## meander1 0.06527705 0.0000 0.0000000 0.0000000 0.0000000 0.0000000
## netcen 0.06898689 0.0000 0.0000000 1.0000000 0.0000000 0.0000000
## updist 0.04496936 0.0000 0.0000000 0.0000000 0.0000000 0.0000000
## BF NA 1.0000 0.6511099 0.3462387 0.2957714 0.2212089
## PostProbs NA 0.6912 0.0450000 0.0239000 0.0204000 0.0153000
## R2 NA 0.0000 0.0906000 0.0565000 0.0478000 0.0315000
## dim NA 1.0000 2.0000000 2.0000000 2.0000000 2.0000000
## logmarg NA 0.0000 -0.4290769 -1.0606268 -1.2181684 -1.5086476image(bas.model, rotate=F)
coef.model <- coef(bas.model)
#abs(coef.model$postmean)-2*coef.model$postsd > 0
plot(confint(coef.model, parm = 2:11))## Warning in arrows(x[not.deg], ci[not.deg, 1], x[not.deg], ci[not.deg, 2], :
## zero-length arrow is of indeterminate angle and so skipped
## NULLconfint <- confint(coef.model, parm = 2:11)
write.table(confint, 'condition.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(3:12),1] <- pip[2:11,1]*sign(coef.model$postmean[2:11])
#coef.model$postmean[2:11]6.2 Variation in host length
bas.model <- bas.lm(avlength ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av +
pool_riffle + meander + netcen + updist,
data=environment2, prior="JZS")
plot(bas.model)
summary(bas.model)## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.0000000 1.0000000 1.00000000 1.0000000 1.000000 1.0000000
## T_av 0.1914598 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## O2_sat_av 0.1407247 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## Con_av 0.4373524 0.0000000 0.00000000 1.0000000 1.000000 0.0000000
## COD_av 0.1445232 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## NH4._av 0.2026061 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## Nt_av 0.1889236 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## pool_riffle1 0.2259504 0.0000000 0.00000000 0.0000000 0.000000 1.0000000
## meander1 0.3217534 0.0000000 0.00000000 0.0000000 1.000000 0.0000000
## netcen 0.6121912 1.0000000 0.00000000 0.0000000 0.000000 1.0000000
## updist 0.1480063 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## BF NA 0.8278694 0.07277547 0.2894406 1.000000 0.6785609
## PostProbs NA 0.1582000 0.13910000 0.0553000 0.042500 0.0288000
## R2 NA 0.2306000 0.00000000 0.1819000 0.314300 0.2982000
## dim NA 2.0000000 1.00000000 2.0000000 3.000000 3.0000000
## logmarg NA 2.4314765 0.00000000 1.3805712 2.620376 2.2325953image(bas.model, rotate=F)
coef.model <- coef(bas.model)
#abs(coef.model$postmean)-2*coef.model$postsd > 0
plot(confint(coef.model, parm = 2:11))## Warning in arrows(x[not.deg], ci[not.deg, 1], x[not.deg], ci[not.deg, 2], :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x[not.deg], ci[not.deg, 1], x[not.deg], ci[not.deg, 2], :
## zero-length arrow is of indeterminate angle and so skipped
## NULLconfint <- confint(coef.model, parm = 2:11)
write.table(confint, 'length.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(3:12),1] <- pip[2:11,1]*sign(coef.model$postmean[2:11])
#coef.model$postmean[2:11]
# Prediction plot
newdata = as.data.frame(cbind(rep(mean(environment2$T_av), 37),
rep(mean(environment2$O2_sat_av), 37),
rep(mean(environment2$Con_av), 37),
rep(mean(environment2$COD_av), 37),
rep(mean(environment2$NH4._av), 37),
rep(mean(environment2$Nt_av), 37),
rep(1, 37),
rep(1, 37),
rep(mean(netcen), 37),
rep(mean(updist), 37)))
colnames(newdata) <- c("T_av", "O2_sat_av", "Con_av", "COD_av", "NH4._av", "Nt_av", "pool_riffle", "meander", "netcen", "updist")
newdata[,"pool_riffle"] <- as.factor(newdata[,"pool_riffle"]); newdata[,"meander"] <- as.factor(newdata[,"meander"])
newdata1 <- newdata; newdata1[,"netcen"] <- netcen
BMA_avlength_netcen <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3j <- ggplot(environment2, aes(netcen, BMA_avlength_netcen$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_avlength_netcen$fit-BMA_avlength_netcen$se.bma.fit), ymax = (BMA_avlength_netcen$fit+BMA_avlength_netcen$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=netcen, y=avlength)) +
labs(x=expression("Network peripherality [m]"), y=expression("Average host length [mm]")) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "j)")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.figure3j
6.3 Variation in Gyrodactylus infection
6.3.1 Mean abundance
bas.model <- bas.lm(avab$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.0000000 1.000000 1.00000000 1.0000000 1.0000000 1.0000000
## avlength 0.5500651 1.000000 0.00000000 0.0000000 0.0000000 1.0000000
## avcondition 0.3049274 0.000000 0.00000000 0.0000000 0.0000000 1.0000000
## T_av 0.1620852 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## O2_sat_av 0.2017699 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## Con_av 0.2621665 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## COD_av 0.8086300 1.000000 0.00000000 1.0000000 1.0000000 1.0000000
## NH4._av 0.2417809 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## Nt_av 0.8618345 1.000000 1.00000000 1.0000000 1.0000000 1.0000000
## pool_riffle1 0.1760484 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## meander1 0.1775842 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## netcen 0.3006370 0.000000 0.00000000 1.0000000 0.0000000 0.0000000
## updist 0.1648814 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## BF NA 1.000000 0.03902505 0.5452757 0.1537678 0.8271829
## PostProbs NA 0.072100 0.05160000 0.0393000 0.0369000 0.0265000
## R2 NA 0.524000 0.28790000 0.5059000 0.4100000 0.5631000
## dim NA 4.000000 2.00000000 4.0000000 3.0000000 5.0000000
## logmarg NA 6.997337 3.75378552 6.3908733 5.1250253 6.8076077image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSEconfint(coef.model)## 2.5% 97.5% beta
## Intercept 9.306741e-02 2.009914e-01 1.489134e-01
## avlength -2.516824e-02 0.000000e+00 -8.122675e-03
## avcondition -1.008942e+00 4.722253e-04 -1.500326e-01
## T_av -1.933113e-02 1.936668e-02 -1.493132e-04
## O2_sat_av -1.010761e-03 5.405876e-03 3.794809e-04
## Con_av -7.801740e-05 5.075870e-04 5.726233e-05
## COD_av 0.000000e+00 1.295780e-02 6.652983e-03
## NH4._av -6.043975e-02 7.670275e-02 -1.888173e-04
## Nt_av 0.000000e+00 5.253835e-02 2.930918e-02
## pool_riffle1 -3.248206e-02 1.064209e-01 3.966550e-03
## meander1 -9.521821e-02 4.524471e-02 -4.806259e-03
## netcen -4.356805e-07 9.961756e-06 1.332021e-06
## updist -1.428611e-06 1.623940e-06 4.618571e-08
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"plot(confint(coef.model, parm = 2:11))
## NULLconfint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])6.3.1.1 Prediction plot for marginal effect of host condition on average Gyrodactylus infection
# Prediction plot
newdata = as.data.frame(cbind(rep(mean(avlength), 37),
rep(mean(avcondition), 37),
rep(mean(environment2$T_av), 37),
rep(mean(environment2$O2_sat_av), 37),
rep(mean(environment2$Con_av), 37),
rep(mean(environment2$COD_av), 37),
rep(mean(environment2$NH4._av), 37),
rep(mean(environment2$Nt_av), 37),
rep(1, 37),
rep(1, 37),
rep(mean(netcen), 37),
rep(mean(updist), 37)))
colnames(newdata) <- c("avlength", "avcondition", "T_av", "O2_sat_av", "Con_av", "COD_av", "NH4._av", "Nt_av", "pool_riffle", "meander", "netcen", "updist")
newdata[,"pool_riffle"] <- as.factor(newdata[,"pool_riffle"]); newdata[,"meander"] <- as.factor(newdata[,"meander"])
newdata1 <- newdata; newdata1[,"avcondition"] <- avcondition
BMA_Gyr_avcond <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3k = ggplot(environment2, aes(avcondition, BMA_Gyr_avcond$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_Gyr_avcond$fit-BMA_Gyr_avcond$se.bma.fit), ymax = (BMA_Gyr_avcond$fit+BMA_Gyr_avcond$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=avcondition, y=avab$Gyr)) +
labs(x=expression("Average host condition"), y=expression(paste(italic("Gyrodactylus"), " spp. - abundance"))) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "k)")
figure3k
6.3.1.2 Prediction plot for marginal effect of COD on average Gyrodactylus infection
newdata1 <- newdata; newdata1[,"COD_av"] <- environment2$COD_av
BMA_Gyr_COD_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3a = ggplot(environment2, aes(COD_av, BMA_Gyr_COD_av$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_Gyr_COD_av$fit-BMA_Gyr_COD_av$se.bma.fit), ymax = (BMA_Gyr_COD_av$fit+BMA_Gyr_COD_av$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=COD_av, y=avab$Gyr)) +
labs(x=expression("Chemical Oxygen Demand [mg/L]"), y=expression(paste(italic("Gyrodactylus"), " spp. - abundance"))) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "a)")
figure3a
6.3.2.3 Prediction plot for marginal effect of total nitrogen on average Gyrodactylus infection
newdata1 <- newdata; newdata1[,"Nt_av"] <- environment2$Nt_av
BMA_Gyr_Nt_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3e = ggplot(environment2, aes(Nt_av, BMA_Gyr_Nt_av$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_Gyr_Nt_av$fit-BMA_Gyr_Nt_av$se.bma.fit), ymax = (BMA_Gyr_Nt_av$fit+BMA_Gyr_Nt_av$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=Nt_av, y=avab$Gyr)) +
labs(x=expression("Nt [mg/L]"), y=expression(paste(italic("Gyrodactylus"), " spp. - abundance"))) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "e)")
figure3e
6.3.2 Median infection intensity
bas.model <- bas.lm(medin$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.00000000 1.0000 1.0000000 1.000000 1.0000000 1.0000000
## avlength 0.02065469 0.0000 0.0000000 0.000000 0.0000000 0.0000000
## avcondition 0.04286905 0.0000 1.0000000 0.000000 0.0000000 0.0000000
## T_av 0.02991678 0.0000 0.0000000 0.000000 0.0000000 0.0000000
## O2_sat_av 0.02400544 0.0000 0.0000000 0.000000 0.0000000 0.0000000
## Con_av 0.02143314 0.0000 0.0000000 0.000000 0.0000000 0.0000000
## COD_av 0.03655302 0.0000 0.0000000 0.000000 1.0000000 0.0000000
## NH4._av 0.03250611 0.0000 0.0000000 0.000000 0.0000000 1.0000000
## Nt_av 0.04088810 0.0000 0.0000000 1.000000 0.0000000 0.0000000
## pool_riffle1 0.02307275 0.0000 0.0000000 0.000000 0.0000000 0.0000000
## meander1 0.02038443 0.0000 0.0000000 0.000000 0.0000000 0.0000000
## netcen 0.02173649 0.0000 0.0000000 0.000000 0.0000000 0.0000000
## updist 0.02738794 0.0000 0.0000000 0.000000 0.0000000 0.0000000
## BF NA 1.0000 0.3176354 0.302852 0.2488803 0.2411373
## PostProbs NA 0.7775 0.0206000 0.019600 0.0161000 0.0156000
## R2 NA 0.0000 0.0517000 0.049100 0.0381000 0.0363000
## dim NA 1.0000 2.0000000 2.000000 2.0000000 2.0000000
## logmarg NA 0.0000 -1.1468512 -1.194511 -1.3907833 -1.4223890image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSEconfint(coef.model)## 2.5% 97.5% beta
## Intercept 1.607225 3.263557 2.405405e+00
## avlength 0.000000 0.000000 -5.938396e-04
## avcondition 0.000000 0.000000 -2.939823e-01
## T_av 0.000000 0.000000 -7.368692e-03
## O2_sat_av 0.000000 0.000000 -4.581964e-04
## Con_av 0.000000 0.000000 -1.594570e-05
## COD_av 0.000000 0.000000 1.609035e-03
## NH4._av 0.000000 0.000000 1.027436e-02
## Nt_av 0.000000 0.000000 7.570653e-03
## pool_riffle1 0.000000 0.000000 1.051497e-02
## meander1 0.000000 0.000000 2.511911e-03
## netcen 0.000000 0.000000 3.990331e-07
## updist 0.000000 0.000000 4.560731e-07
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"plot(confint(coef.model, parm = 2:11))
## NULLconfint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])6.3.2 Prevalence
bas.model <- bas.lm(prev$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.00000000 1.000000 1.00000000 1.0000000 1.0000000 1.00000000
## avlength 0.16685753 0.000000 0.00000000 0.0000000 1.0000000 0.00000000
## avcondition 0.07728991 0.000000 0.00000000 0.0000000 0.0000000 0.00000000
## T_av 0.07242807 0.000000 0.00000000 0.0000000 0.0000000 0.00000000
## O2_sat_av 0.09667222 0.000000 0.00000000 0.0000000 0.0000000 0.00000000
## Con_av 0.58076129 1.000000 0.00000000 1.0000000 0.0000000 0.00000000
## COD_av 0.17268847 0.000000 0.00000000 1.0000000 0.0000000 0.00000000
## NH4._av 0.09134534 0.000000 0.00000000 0.0000000 0.0000000 0.00000000
## Nt_av 0.11056039 0.000000 0.00000000 0.0000000 0.0000000 0.00000000
## pool_riffle1 0.13922799 0.000000 0.00000000 0.0000000 0.0000000 0.00000000
## meander1 0.09905733 0.000000 0.00000000 0.0000000 0.0000000 0.00000000
## netcen 0.11864979 0.000000 0.00000000 0.0000000 0.0000000 1.00000000
## updist 0.06837623 0.000000 0.00000000 0.0000000 0.0000000 0.00000000
## BF NA 1.000000 0.08283414 0.8857497 0.1246735 0.07409002
## PostProbs NA 0.228800 0.22750000 0.0369000 0.0285000 0.01700000
## R2 NA 0.233300 0.00000000 0.3039000 0.1341000 0.10730000
## dim NA 2.000000 1.00000000 3.0000000 2.0000000 2.00000000
## logmarg NA 2.490915 0.00000000 2.3695941 0.4088580 -0.11155941image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSEconfint(coef.model)## 2.5% 97.5% beta
## Intercept 2.547614e-01 3.964116e-01 3.224082e-01
## avlength -2.300336e-02 0.000000e+00 -2.431148e-03
## avcondition -3.425493e-01 9.865432e-02 -1.986363e-02
## T_av -1.891814e-02 2.173153e-04 -7.123130e-04
## O2_sat_av -4.068781e-03 1.262431e-04 -2.662986e-04
## Con_av 0.000000e+00 8.064140e-04 3.112530e-04
## COD_av 0.000000e+00 8.036852e-03 9.376253e-04
## NH4._av -3.593614e-03 4.398906e-02 1.064409e-03
## Nt_av -1.789721e-04 2.438631e-02 1.612509e-03
## pool_riffle1 0.000000e+00 1.577918e-01 1.465705e-02
## meander1 -1.043718e-01 4.991691e-04 -6.966751e-03
## netcen -2.974468e-09 8.169806e-06 5.765326e-07
## updist -8.239812e-07 1.990176e-07 -6.922933e-09
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"plot(confint(coef.model, parm = 2:11))
## NULLconfint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])6.3.2.3 Prediction plot for marginal effect of conductivity on average Gyrodactylus infection
newdata1 <- newdata; newdata1[,"Con_av"] <- environment2$Con_av
BMA_Gyr_Con_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3i = ggplot(environment2, aes(Con_av, BMA_Gyr_Con_av$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_Gyr_Con_av$fit-BMA_Gyr_Con_av$se.bma.fit), ymax = (BMA_Gyr_Con_av$fit+BMA_Gyr_Con_av$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=Con_av, y=prev$Gyr)) +
labs(x=expression(paste("Conductivity [", mu, "S/cm]")), y=expression(paste(italic("Gyrodactylus"), " spp. - prevalence"))) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "i)")
figure3i
6.4 Variation in Trichodina infection
6.4.1 Mean abundance
bas.model <- bas.lm(avab$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.00000000 1.0000000 1.000000 1.000000 1.0000000 1.0000000
## avlength 0.12444042 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## avcondition 0.09559152 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## T_av 0.08500185 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## O2_sat_av 0.09751169 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## Con_av 0.42180029 0.0000000 1.000000 0.000000 0.0000000 0.0000000
## COD_av 0.46724539 0.0000000 1.000000 0.000000 0.0000000 1.0000000
## NH4._av 0.21465122 0.0000000 0.000000 0.000000 1.0000000 0.0000000
## Nt_av 0.34576751 0.0000000 0.000000 1.000000 0.0000000 1.0000000
## pool_riffle1 0.13743699 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## meander1 0.08968689 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## netcen 0.10313787 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## updist 0.08604777 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## BF NA 0.0239404 1.000000 0.170643 0.1105942 0.2711811
## PostProbs NA 0.1650000 0.104400 0.098000 0.0635000 0.0283000
## R2 NA 0.0000000 0.358500 0.209300 0.1890000 0.3063000
## dim NA 1.0000000 3.000000 2.000000 2.0000000 3.0000000
## logmarg NA 0.0000000 3.732188 1.964006 1.5303004 2.4272192image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSEconfint(coef.model)## 2.5% 97.5% beta
## Intercept 1.232338e-01 2.714579e-01 2.002413e-01
## avlength -1.705124e-02 3.780264e-04 -1.114390e-03
## avcondition -5.882325e-01 5.136492e-02 -2.394139e-02
## T_av -1.537747e-02 1.020332e-02 -1.479738e-04
## O2_sat_av -1.352662e-03 2.779387e-03 1.132458e-04
## Con_av 0.000000e+00 7.626598e-04 2.109076e-04
## COD_av -1.738195e-05 1.373697e-02 4.203136e-03
## NH4._av -8.071803e-04 1.061136e-01 1.096330e-02
## Nt_av 0.000000e+00 4.926715e-02 1.044647e-02
## pool_riffle1 0.000000e+00 1.494898e-01 1.212839e-02
## meander1 -7.138422e-02 2.663792e-02 -1.678368e-03
## netcen -6.434350e-06 4.469308e-07 -2.846749e-07
## updist -1.450572e-06 7.187041e-07 -1.915234e-08
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"plot(confint(coef.model, parm = 2:11))
## NULLconfint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])6.4.1 Median infection intensity
bas.model <- bas.lm(medin$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.0000000 1.000000 1.00000000 1.0000000 1.0000000 1.0000000
## avlength 0.2881774 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## avcondition 0.2224163 0.000000 0.00000000 0.0000000 1.0000000 0.0000000
## T_av 0.1574615 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## O2_sat_av 0.1621652 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## Con_av 0.7912882 1.000000 0.00000000 1.0000000 1.0000000 1.0000000
## COD_av 0.7219333 1.000000 0.00000000 0.0000000 1.0000000 1.0000000
## NH4._av 0.3248706 0.000000 1.00000000 1.0000000 0.0000000 0.0000000
## Nt_av 0.3218716 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## pool_riffle1 0.2081402 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## meander1 0.2655250 0.000000 0.00000000 0.0000000 0.0000000 1.0000000
## netcen 0.2107764 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## updist 0.1536937 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## BF NA 1.000000 0.04957628 0.2279222 0.5092549 0.4837322
## PostProbs NA 0.148800 0.04060000 0.0339000 0.0227000 0.0216000
## R2 NA 0.449900 0.26810000 0.3987000 0.4816000 0.4799000
## dim NA 3.000000 2.00000000 3.0000000 4.0000000 4.0000000
## logmarg NA 6.288599 3.28435583 4.8098478 5.6137920 5.5623749image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSEconfint(coef.model)## 2.5% 97.5% beta
## Intercept 1.675518e+01 3.441505e+01 2.547297e+01
## avlength -3.471106e+00 1.978629e-01 -5.158998e-01
## avcondition -1.452798e+02 8.422586e-02 -1.375058e+01
## T_av -1.669965e+00 4.431680e+00 2.065105e-01
## O2_sat_av -4.995746e-01 4.824268e-01 -1.498026e-02
## Con_av 0.000000e+00 1.283962e-01 6.585727e-02
## COD_av 0.000000e+00 2.001430e+00 9.195418e-01
## NH4._av -2.433043e+00 1.595274e+01 2.222211e+00
## Nt_av -3.006566e-01 6.163537e+00 1.012166e+00
## pool_riffle1 -2.712539e+00 2.385458e+01 2.247962e+00
## meander1 -2.677652e+01 2.623704e-01 -3.523477e+00
## netcen -1.459503e-03 1.344122e-04 -1.366046e-04
## updist -2.059586e-04 2.662744e-04 5.653370e-06
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"plot(confint(coef.model, parm = 2:11))
## NULLconfint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])6.3.2.3 Prediction plot for marginal effect of conductivity on median infection intensity with Trichodina
newdata1 <- newdata; newdata1[,"Con_av"] <- environment2$Con_av
BMA_Tri_Con_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3h = ggplot(environment2, aes(Con_av, BMA_Tri_Con_av$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_Tri_Con_av$fit-BMA_Tri_Con_av$se.bma.fit), ymax = (BMA_Tri_Con_av$fit+BMA_Tri_Con_av$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=Con_av, y=medin$Tri)) +
labs(x=expression(paste("Conductivity [", mu, "S/cm]")), y=expression(paste(italic("Trichodina"), " sp. - infection intensity"))) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "h)")
figure3h
6.3.2.3 Prediction plot for marginal effect of COD on median infection intensity with Trichodina
newdata1 <- newdata; newdata1[,"COD_av"] <- environment2$COD_av
BMA_Tri_COD_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3d = ggplot(environment2, aes(COD_av, BMA_Tri_COD_av$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_Tri_COD_av$fit-BMA_Tri_COD_av$se.bma.fit), ymax = (BMA_Tri_COD_av$fit+BMA_Tri_COD_av$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=COD_av, y=medin$Tri)) +
labs(x=expression(paste("Chemical Oxygen Demand [mg/L]")), y=expression(paste(italic("Trichodina"), " sp. - infection intensity"))) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "d)")
figure3d
6.4.1 Prevalence
bas.model <- bas.lm(prev$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000
## avlength 0.04005815 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
## avcondition 0.03389439 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
## T_av 0.03654853 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
## O2_sat_av 0.03749123 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
## Con_av 0.20573471 0.0000000 1.0000000 0.0000000 1.0000000 0.000000
## COD_av 0.04629842 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
## NH4._av 0.06690739 0.0000000 0.0000000 1.0000000 0.0000000 0.000000
## Nt_av 0.04604656 0.0000000 0.0000000 0.0000000 0.0000000 1.000000
## pool_riffle1 0.03411173 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
## meander1 0.03657952 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
## netcen 0.07467299 0.0000000 0.0000000 0.0000000 1.0000000 0.000000
## updist 0.07352673 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
## BF NA 0.5616705 0.7274201 0.2732671 1.0000000 0.155205
## PostProbs NA 0.6354000 0.0686000 0.0258000 0.0171000 0.014600
## R2 NA 0.0000000 0.1264000 0.0750000 0.2249000 0.044000
## dim NA 1.0000000 2.0000000 2.0000000 3.0000000 2.000000
## logmarg NA 0.0000000 0.2585887 -0.7204656 0.5768398 -1.286169image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSEconfint(coef.model)## 2.5% 97.5% beta
## Intercept 7.820697e-01 0.8851797650 8.340580e-01
## avlength 0.000000e+00 0.0000000000 -1.453558e-04
## avcondition 0.000000e+00 0.0000000000 -2.610917e-03
## T_av 0.000000e+00 0.0000000000 -1.135202e-05
## O2_sat_av 0.000000e+00 0.0000000000 -1.893466e-05
## Con_av 0.000000e+00 0.0004278091 6.570599e-05
## COD_av 0.000000e+00 0.0000000000 9.451636e-05
## NH4._av 0.000000e+00 0.0191693460 1.713706e-03
## Nt_av -1.240461e-04 0.0000000000 1.739937e-04
## pool_riffle1 0.000000e+00 0.0000000000 2.104077e-04
## meander1 0.000000e+00 0.0000000000 -6.652944e-04
## netcen -3.878705e-06 0.0000000000 -3.317353e-07
## updist -1.482392e-06 0.0000000000 -1.309710e-07
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"plot(confint(coef.model, parm = 2:11))
## NULLconfint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])6.5 Variation in Glugea infection
bas.model <- bas.glm(pa$Glu ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av +
NH4._av + Nt_av + SM_av + pool_riffle + meander +
spavar$netcen + spavar$updist, data=environment2, betaprior=g.prior(100), family=binomial)
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)## P(B != 0 | Y) model 1 model 2 model 3 model 4
## Intercept 1.0000000 1.000000 1.0000000 1.000000 1.0000000
## avlength 0.6661011 0.000000 1.0000000 0.000000 0.0000000
## avcondition 0.4341431 0.000000 1.0000000 1.000000 1.0000000
## T_av 0.2906006 0.000000 0.0000000 0.000000 0.0000000
## O2_sat_av 0.4036499 0.000000 0.0000000 0.000000 0.0000000
## Con_av 0.9885742 1.000000 1.0000000 1.000000 1.0000000
## COD_av 0.5252563 1.000000 1.0000000 0.000000 1.0000000
## NH4._av 0.8087402 1.000000 1.0000000 1.000000 1.0000000
## Nt_av 0.5433594 1.000000 0.0000000 1.000000 1.0000000
## SM_av 0.4432983 1.000000 0.0000000 1.000000 0.0000000
## pool_riffle1 0.4388672 0.000000 0.0000000 0.000000 0.0000000
## meander1 0.9899902 1.000000 1.0000000 1.000000 1.0000000
## spavar$netcen 0.4765747 0.000000 0.0000000 1.000000 0.0000000
## spavar$updist 0.8123047 1.000000 1.0000000 0.000000 1.0000000
## BF NA 1.000000 0.8189871 0.769804 0.7351301
## PostProbs NA 0.046100 0.0397000 0.038800 0.0341000
## R2 NA 1.000000 1.0000000 1.000000 1.0000000
## dim NA 8.000000 8.0000000 8.000000 8.0000000
## logmarg NA -5.808403 -6.0080902 -6.070023 -6.1161110
## model 5
## Intercept 1.0000000
## avlength 1.0000000
## avcondition 0.0000000
## T_av 0.0000000
## O2_sat_av 1.0000000
## Con_av 1.0000000
## COD_av 0.0000000
## NH4._av 0.0000000
## Nt_av 0.0000000
## SM_av 0.0000000
## pool_riffle1 1.0000000
## meander1 1.0000000
## spavar$netcen 1.0000000
## spavar$updist 1.0000000
## BF 0.5830791
## PostProbs 0.0313000
## R2 1.0000000
## dim 8.0000000
## logmarg -6.3478357image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSEconfint(coef.model)## 2.5% 97.5% beta
## Intercept -2.902737e+06 2.926170e+06 -6.339692e+03
## avlength -4.505021e+04 4.199676e+04 9.029383e+01
## avcondition -9.023933e+05 8.924758e+05 1.496115e+01
## T_av -2.790922e+04 3.377942e+04 1.699320e+01
## O2_sat_av -5.691607e+03 6.068530e+03 9.150769e+00
## Con_av -1.608333e+03 1.623631e+03 4.293194e+00
## COD_av -1.527908e+04 1.311479e+04 9.080887e+00
## NH4._av -2.332931e+05 2.150607e+05 -2.537950e+02
## Nt_av -4.507645e+04 4.741178e+04 3.593276e+01
## SM_av -2.247258e+03 2.416569e+03 -1.167050e+00
## pool_riffle1 -1.421694e+05 1.397596e+05 8.653441e+01
## meander1 -4.926730e+05 5.194216e+05 -1.344934e+03
## spavar$netcen -1.058765e+01 9.300020e+00 -8.848165e-03
## spavar$updist -4.812767e+00 5.076095e+00 -7.824641e-03
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"plot(confint(coef.model, parm = 2:11))
## NULLconfint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),9] <- pip[2:13,1]*sign(coef.model$postmean[2:13])6.5 Variation in Contracaecum infection
bas.model <- bas.glm(pa$Con ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av +
NH4._av + Nt_av + SM_av + pool_riffle + meander +
spavar$netcen + spavar$updist, data=environment2, betaprior=g.prior(100), family=binomial)
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)## P(B != 0 | Y) model 1 model 2 model 3 model 4
## Intercept 1.00000000 1.00000 1.0000000 1.0000000 1.0000000
## avlength 0.08469238 0.00000 1.0000000 0.0000000 0.0000000
## avcondition 0.02604980 0.00000 0.0000000 1.0000000 0.0000000
## T_av 0.01040039 0.00000 0.0000000 0.0000000 0.0000000
## O2_sat_av 0.02716064 0.00000 0.0000000 0.0000000 1.0000000
## Con_av 0.01129150 0.00000 0.0000000 0.0000000 0.0000000
## COD_av 0.01397705 0.00000 0.0000000 0.0000000 0.0000000
## NH4._av 0.01912842 0.00000 0.0000000 0.0000000 0.0000000
## Nt_av 0.01209717 0.00000 0.0000000 0.0000000 0.0000000
## SM_av 0.01187744 0.00000 0.0000000 0.0000000 0.0000000
## pool_riffle1 0.01187744 0.00000 0.0000000 0.0000000 0.0000000
## meander1 0.01021729 0.00000 0.0000000 0.0000000 0.0000000
## spavar$netcen 0.01888428 0.00000 0.0000000 0.0000000 0.0000000
## spavar$updist 0.01214600 0.00000 0.0000000 0.0000000 0.0000000
## BF NA 1.00000 0.9507459 0.2581488 0.2583995
## PostProbs NA 0.78470 0.0608000 0.0152000 0.0145000
## R2 NA 0.00000 0.0864000 0.0365000 0.0366000
## dim NA 1.00000 2.0000000 2.0000000 2.0000000
## logmarg NA -25.82594 -25.8764468 -27.1801577 -27.1791867
## model 5
## Intercept 1.00000
## avlength 0.00000
## avcondition 0.00000
## T_av 0.00000
## O2_sat_av 0.00000
## Con_av 0.00000
## COD_av 0.00000
## NH4._av 1.00000
## Nt_av 0.00000
## SM_av 0.00000
## pool_riffle1 0.00000
## meander1 0.00000
## spavar$netcen 0.00000
## spavar$updist 0.00000
## BF 0.19791
## PostProbs 0.01250
## R2 0.02640
## dim 2.00000
## logmarg -27.44588image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSEconfint(coef.model)## 2.5% 97.5% beta
## Intercept -9.938373 2.1628712 -5.683824e-01
## avlength 0.000000 0.1692877 1.573802e-02
## avcondition 0.000000 0.0000000 -1.691605e-01
## T_av 0.000000 0.0000000 7.921083e-04
## O2_sat_av 0.000000 0.0000000 1.190308e-03
## Con_av 0.000000 0.0000000 5.101715e-06
## COD_av 0.000000 0.0000000 -3.712753e-04
## NH4._av 0.000000 0.0000000 -5.964150e-03
## Nt_av 0.000000 0.0000000 -9.152172e-04
## SM_av 0.000000 0.0000000 1.162131e-04
## pool_riffle1 0.000000 0.0000000 -5.259321e-03
## meander1 0.000000 0.0000000 2.907684e-03
## spavar$netcen 0.000000 0.0000000 -8.702274e-07
## spavar$updist 0.000000 0.0000000 -1.184043e-07
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"plot(confint(coef.model, parm = 2:11))
## NULLconfint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),10] <- pip[2:13,1]*sign(coef.model$postmean[2:13])Variation in Anguillicoloides infection
bas.model <- bas.glm(pa$Ang ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av +
NH4._av + Nt_av + SM_av + pool_riffle + meander +
spavar$netcen + spavar$updist, data=environment2, betaprior=g.prior(100), family=binomial)
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)## P(B != 0 | Y) model 1 model 2 model 3 model 4
## Intercept 1.0000000 1.000000 1.0000000 1.0000000 1.0000000
## avlength 0.9313965 1.000000 1.0000000 1.0000000 1.0000000
## avcondition 0.3154663 0.000000 0.0000000 0.0000000 0.0000000
## T_av 0.9406494 1.000000 1.0000000 1.0000000 1.0000000
## O2_sat_av 0.8143555 1.000000 1.0000000 1.0000000 1.0000000
## Con_av 0.5466187 0.000000 1.0000000 1.0000000 0.0000000
## COD_av 0.3676270 0.000000 0.0000000 1.0000000 0.0000000
## NH4._av 0.9558350 1.000000 1.0000000 1.0000000 1.0000000
## Nt_av 0.5304565 0.000000 1.0000000 0.0000000 0.0000000
## SM_av 0.9200195 1.000000 1.0000000 1.0000000 1.0000000
## pool_riffle1 0.3712646 0.000000 0.0000000 0.0000000 1.0000000
## meander1 0.9378540 1.000000 1.0000000 1.0000000 1.0000000
## spavar$netcen 0.8593750 1.000000 1.0000000 1.0000000 1.0000000
## spavar$updist 0.5289795 1.000000 0.0000000 0.0000000 1.0000000
## BF NA 1.000000 0.4728289 0.3794633 0.2858158
## PostProbs NA 0.079300 0.0720000 0.0531000 0.0442000
## R2 NA 1.000000 1.0000000 1.0000000 1.0000000
## dim NA 9.000000 10.0000000 10.0000000 10.0000000
## logmarg NA -9.802499 -10.5515211 -10.7714968 -11.0549070
## model 5
## Intercept 1.0000000
## avlength 1.0000000
## avcondition 1.0000000
## T_av 1.0000000
## O2_sat_av 0.0000000
## Con_av 1.0000000
## COD_av 0.0000000
## NH4._av 1.0000000
## Nt_av 1.0000000
## SM_av 1.0000000
## pool_riffle1 0.0000000
## meander1 1.0000000
## spavar$netcen 1.0000000
## spavar$updist 0.0000000
## BF 0.2037316
## PostProbs 0.0299000
## R2 1.0000000
## dim 10.0000000
## logmarg -11.3934513image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSEconfint(coef.model)## 2.5% 97.5% beta
## Intercept -2.493583e+06 2.220172e+06 -6.314576e+03
## avlength -4.633181e+04 4.189389e+04 1.289938e+02
## avcondition -8.259937e+05 7.290298e+05 -2.425364e+02
## T_av -1.617354e+05 1.522072e+05 5.783392e+02
## O2_sat_av -1.894797e+04 1.783841e+04 -3.694674e+01
## Con_av -6.725601e+02 6.298949e+02 3.253851e-01
## COD_av -7.695858e+03 7.302203e+03 1.422581e+00
## NH4._av -3.522546e+05 3.350507e+05 -1.031216e+03
## Nt_av -4.309152e+04 4.253718e+04 -4.515461e+01
## SM_av -5.880688e+03 5.728770e+03 1.263286e+01
## pool_riffle1 -1.148759e+05 1.313316e+05 -1.031862e+02
## meander1 -3.968940e+05 3.849816e+05 -1.367302e+03
## spavar$netcen -1.525973e+01 1.372595e+01 -2.731076e-02
## spavar$updist -3.007013e+00 2.830525e+00 1.319948e-03
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"plot(confint(coef.model, parm = 2:11))
## NULLconfint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),11] <- pip[2:13,1]*sign(coef.model$postmean[2:13])Variation in Individual Parasitation Index (all parasites)
bas.model <- bas.lm(avPI ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.0000000 1.00000000 1.00000000 1.00000000 1.000000 1.000000
## avlength 0.2261960 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## avcondition 0.1810890 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## T_av 0.4574073 0.00000000 1.00000000 0.00000000 1.000000 0.000000
## O2_sat_av 0.1952767 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## Con_av 0.4598919 0.00000000 0.00000000 0.00000000 0.000000 1.000000
## COD_av 0.5684747 0.00000000 0.00000000 0.00000000 0.000000 1.000000
## NH4._av 0.2748144 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## Nt_av 0.4582802 0.00000000 0.00000000 1.00000000 1.000000 0.000000
## pool_riffle1 0.3054914 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## meander1 0.1985617 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## netcen 0.1910872 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## updist 0.4153657 0.00000000 0.00000000 0.00000000 0.000000 1.000000
## BF NA 0.01969453 0.09251289 0.08737546 0.355023 1.000000
## PostProbs NA 0.10850000 0.04250000 0.04010000 0.029600 0.025000
## R2 NA 0.00000000 0.18980000 0.18710000 0.325400 0.424800
## dim NA 1.00000000 2.00000000 2.00000000 3.000000 4.000000
## logmarg NA 0.00000000 1.54700725 1.48987379 2.891842 3.927415image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSEconfint(coef.model)## 2.5% 97.5% beta
## Intercept 1.139112e+00 1.732523e+00 1.422090e+00
## avlength -8.655588e-02 1.877262e-02 -7.992245e-03
## avcondition -2.976576e+00 1.529684e+00 -9.325811e-02
## T_av -7.038395e-04 3.716692e-01 8.644408e-02
## O2_sat_av -2.345726e-02 1.521722e-02 -7.476618e-04
## Con_av -3.164277e-06 3.173245e-03 8.405554e-04
## COD_av 0.000000e+00 6.148546e-02 2.109188e-02
## NH4._av -5.306465e-01 1.492837e-01 -4.965441e-02
## Nt_av 0.000000e+00 2.234194e-01 5.507750e-02
## pool_riffle1 -3.737605e-02 8.953177e-01 1.389140e-01
## meander1 -5.590964e-01 2.250265e-01 -3.866036e-02
## netcen -2.998202e-05 1.708883e-05 -1.270941e-06
## updist -2.627701e-05 0.000000e+00 -5.952452e-06
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"plot(confint(coef.model, parm = 2:11))
## NULLconfint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),12] <- pip[2:13,1]*sign(coef.model$postmean[2:13])Prediction plot for marginal effect of COD on Individual Parasitation Index
newdata1 <- newdata; newdata1[,"COD_av"] <- environment2$COD_av
BMA_IPI_COD_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3b = ggplot(environment2, aes(COD_av, BMA_IPI_COD_av$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_IPI_COD_av$fit-BMA_IPI_COD_av$se.bma.fit), ymax = (BMA_IPI_COD_av$fit+BMA_IPI_COD_av$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=COD_av, y=avPI)) +
labs(x=expression("Chemical Oxygen Demand [mg/L]"), y=expression("I"[PI]*" - all parasites")) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "b)")
figure3b
Variation in Individual Parasitation Index (only ectoparasites)
bas.model <- bas.lm(avPI_ecto ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.0000000 1.0000000 1.000000000 1.0000000 1.000000 1.000000
## avlength 0.4669536 0.0000000 1.000000000 0.0000000 0.000000 1.000000
## avcondition 0.3758074 0.0000000 1.000000000 0.0000000 0.000000 0.000000
## T_av 0.2753561 0.0000000 1.000000000 0.0000000 0.000000 0.000000
## O2_sat_av 0.2919614 0.0000000 1.000000000 0.0000000 0.000000 0.000000
## Con_av 0.6617062 0.0000000 1.000000000 1.0000000 0.000000 0.000000
## COD_av 0.8437205 0.0000000 1.000000000 1.0000000 1.000000 1.000000
## NH4._av 0.4163705 0.0000000 1.000000000 0.0000000 0.000000 0.000000
## Nt_av 0.8302866 1.0000000 1.000000000 0.0000000 1.000000 1.000000
## pool_riffle1 0.4025108 0.0000000 1.000000000 0.0000000 0.000000 0.000000
## meander1 0.5031625 0.0000000 1.000000000 0.0000000 0.000000 0.000000
## netcen 0.3010717 0.0000000 1.000000000 0.0000000 0.000000 0.000000
## updist 0.3285744 0.0000000 1.000000000 0.0000000 0.000000 0.000000
## BF NA 0.1318791 0.008313608 0.4282678 0.319205 1.000000
## PostProbs NA 0.0510000 0.038600000 0.0301000 0.022500 0.021100
## R2 NA 0.3007000 0.682800000 0.4141000 0.403700 0.496600
## dim NA 2.0000000 13.000000000 3.0000000 3.000000 4.000000
## logmarg NA 4.0635523 1.299560343 5.2414154 4.947500 6.089422image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSEconfint(coef.model)## 2.5% 97.5% beta
## Intercept 6.885848e-01 1.160730e+00 9.199577e-01
## avlength -1.053581e-01 1.014527e-02 -2.357999e-02
## avcondition -4.323979e+00 6.216264e-01 -7.055549e-01
## T_av -1.006276e-01 1.379698e-01 5.818763e-03
## O2_sat_av -1.303563e-02 2.535981e-02 1.122707e-03
## Con_av -2.688244e-06 3.154139e-03 1.203349e-03
## COD_av 0.000000e+00 6.418310e-02 3.405818e-02
## NH4._av -4.708824e-01 1.076874e-01 -8.034086e-02
## Nt_av 0.000000e+00 2.533901e-01 1.257830e-01
## pool_riffle1 -8.969849e-02 7.899060e-01 1.490513e-01
## meander1 -9.360961e-01 1.708101e-02 -2.358742e-01
## netcen -3.378862e-05 2.125678e-05 -2.447510e-06
## updist -1.563006e-05 4.018463e-06 -1.836166e-06
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"plot(confint(coef.model, parm = 2:11))
## NULLconfint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),13] <- pip[2:13,1]*sign(coef.model$postmean[2:13])Prediction plot for marginal effect of COD on Individual Parasitation Index (only ectoparasites)
newdata1 <- newdata; newdata1[,"COD_av"] <- environment2$COD_av
BMA_IPIecto_COD_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3c = ggplot(environment2, aes(COD_av, BMA_IPIecto_COD_av$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_IPIecto_COD_av$fit-BMA_IPIecto_COD_av$se.bma.fit), ymax = (BMA_IPIecto_COD_av$fit+BMA_IPIecto_COD_av$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=COD_av, y=avPI_ecto)) +
labs(x=expression("Chemical Oxygen Demand [mg/L]"), y=expression("I"[PI]*" - ectoparasites")) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "c)")
figure3c
Prediction plot for marginal effect of total nitrogen on Individual Parasitation Index (only ectoparasites)
newdata1 <- newdata; newdata1[,"Nt_av"] <- environment2$Nt_av
BMA_IPIecto_Nt_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3f = ggplot(environment2, aes(Nt_av, BMA_IPIecto_Nt_av$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_IPIecto_Nt_av$fit-BMA_IPIecto_Nt_av$se.bma.fit), ymax = (BMA_IPIecto_Nt_av$fit+BMA_IPIecto_Nt_av$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=Nt_av, y=avPI_ecto)) +
labs(x=expression("Nt [mg/L]"), y=expression("I"[PI]*" - ectoparasites")) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "f)")
figure3f
Prediction plot for marginal effect of conductivity on Individual Parasitation Index (only ectoparasites)
newdata1 <- newdata; newdata1[,"Con_av"] <- environment2$Con_av
BMA_IPIecto_Con_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3g = ggplot(environment2, aes(Con_av, BMA_IPIecto_Con_av$fit)) +
theme_bw() +
geom_line(color="red", size=1) +
geom_ribbon(aes(ymin = (BMA_IPIecto_Con_av$fit-BMA_IPIecto_Con_av$se.bma.fit), ymax = (BMA_IPIecto_Con_av$fit+BMA_IPIecto_Con_av$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=Con_av, y=avPI)) +
labs(x=expression(paste("Conductivity [", mu, "S/cm]")), y=expression("I"[PI]*" - all parasites")) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "g)")
figure3g
Prediction plot for marginal effect of meanders on Individual Parasitation Index (only ectoparasites)
newdata1 <- newdata; newdata1[,"meander"] <- environment2$meander
BMA_IPIecto_meander <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3l = ggplot(environment2, aes(meander, BMA_IPIecto_meander$fit)) +
theme_bw() +
#geom_line(color="red", size=1) +
#geom_ribbon(aes(ymin = (BMA$fit-BMA$se.bma.fit), ymax = (BMA$fit+BMA$se.bma.fit)), alpha = .1) +
geom_point(data = environment2, aes(x=meander, y=avPI)) +
geom_boxplot(aes(lower = (BMA_IPIecto_meander$fit-BMA_IPIecto_meander$se.bma.fit), middle = BMA_IPIecto_meander$fit, upper = (BMA_IPIecto_meander$fit+BMA_IPIecto_meander$se.bma.fit))) +
labs(x=expression("Meander"), y=expression("I"[PI]*" - ectoparasites")) +
theme(axis.title.x = element_text(size=12),
axis.title.y = element_text(size=12)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(subtitle = "l)")
figure3l
Variation in Individual Parasitation Index (only endoparasites)
bas.model <- bas.lm(avPI_endo ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)
summary(bas.model)## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.00000000 1.0000000 1.000000 1.00000000 1.0000000 1.0000000
## avlength 0.06442140 0.0000000 0.000000 0.00000000 0.0000000 0.0000000
## avcondition 0.12246410 0.0000000 0.000000 1.00000000 0.0000000 0.0000000
## T_av 0.30397004 0.0000000 1.000000 0.00000000 0.0000000 1.0000000
## O2_sat_av 0.05178919 0.0000000 0.000000 0.00000000 0.0000000 0.0000000
## Con_av 0.05158494 0.0000000 0.000000 0.00000000 0.0000000 0.0000000
## COD_av 0.04593088 0.0000000 0.000000 0.00000000 0.0000000 0.0000000
## NH4._av 0.08528364 0.0000000 0.000000 0.00000000 0.0000000 1.0000000
## Nt_av 0.06477759 0.0000000 0.000000 0.00000000 0.0000000 0.0000000
## pool_riffle1 0.04504579 0.0000000 0.000000 0.00000000 0.0000000 0.0000000
## meander1 0.05899866 0.0000000 0.000000 0.00000000 0.0000000 0.0000000
## netcen 0.05414430 0.0000000 0.000000 0.00000000 0.0000000 0.0000000
## updist 0.10936917 0.0000000 0.000000 0.00000000 1.0000000 0.0000000
## BF NA 0.3484238 1.000000 0.33059983 0.2877395 0.6134470
## PostProbs NA 0.4921000 0.117700 0.03890000 0.0339000 0.0131000
## R2 NA 0.0000000 0.166100 0.11040000 0.1032000 0.2244000
## dim NA 1.0000000 2.000000 2.00000000 2.0000000 3.0000000
## logmarg NA 0.0000000 1.054336 -0.05251095 -0.1913639 0.5656742image(bas.model, rotate=F)
coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSEconfint(coef.model)## 2.5% 97.5% beta
## Intercept 3.062561e-01 6.936441e-01 4.947994e-01
## avlength -1.350725e-04 1.781226e-02 1.398963e-03
## avcondition 0.000000e+00 2.576839e+00 2.446911e-01
## T_av 0.000000e+00 2.160719e-01 4.520999e-02
## O2_sat_av -3.667475e-04 1.264879e-04 -2.446396e-04
## Con_av 0.000000e+00 0.000000e+00 -1.385311e-05
## COD_av 0.000000e+00 0.000000e+00 -1.072193e-04
## NH4._av -1.139164e-01 0.000000e+00 -9.322620e-03
## Nt_av -2.667990e-02 2.022951e-04 -2.149910e-03
## pool_riffle1 0.000000e+00 0.000000e+00 2.267863e-03
## meander1 0.000000e+00 8.606124e-02 1.005422e-02
## netcen -1.545415e-07 1.173066e-06 -3.607730e-07
## updist -8.998181e-06 0.000000e+00 -8.091451e-07
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"plot(confint(coef.model, parm = 2:11))## Warning in arrows(x[not.deg], ci[not.deg, 1], x[not.deg], ci[not.deg, 2], :
## zero-length arrow is of indeterminate angle and so skipped
## NULLconfint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),14] <- pip[2:13,1]*sign(coef.model$postmean[2:13])Figure 3
library(gridExtra)
pdf("Figure3.pdf", width = 10.8, height = 14.4)
grid.arrange(figure3a, figure3b, figure3c, figure3d, figure3e, figure3f, figure3g, figure3h, figure3i, figure3j, figure3k, figure3l)
dev.off()## png
## 27. BORAL analysis
Model-based analysis of multivariate abundance data using Bayesian Markov chain Monte Carlo methods for parameter estimation
7.1 BORAL analysis for average abundances of parasites
data$Site <- as.factor(data$site)
levels(data$site) <- levels(as.factor(environment2$Site))
data_m <- merge(data, environment2, by = "Site")
data_all <- na.omit(data_m)
names(data_all)## [1] "Site" "site" "fish"
## [4] "parspeciesrichness" "div_shannon" "div_simpson"
## [7] "temp" "pH" "conductivity"
## [10] "nitrogen" "phosphorus" "oxygen"
## [13] "netcen.x" "updist.x" "updist2"
## [16] "updist3" "fishspeciesrichness" "weight"
## [19] "weigh..g." "length" "SMI"
## [22] "Sex" "Gyr" "Tri"
## [25] "Glu" "ecto_screener" "Con"
## [28] "CysL" "Pro" "Aca"
## [31] "Cam" "Ang" "CysI"
## [34] "endo_screener" "PI" "PI_ecto"
## [37] "PI_endo" "T_av" "O2_sat_av"
## [40] "Con_av" "COD_av" "NH4._av"
## [43] "Nt_av" "SM_av" "pool_riffle"
## [46] "meander" "updist.y" "netcen.y"avcondition <- aggregate(data$SMI, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2]
avlength <- aggregate(data$length, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2]
y <- round(cbind(avab$Gyr, avab$Tri, avab$Glu, avab$Con, avab$Ang))
X <- cbind(avcondition,
avlength,
environment2$T_av,
environment2$O2_sat_av,
environment2$Con_av,
environment2$COD_av,
environment2$NH4._av,
environment2$Nt_av,
environment2$netcen,
environment2$updist,
as.numeric(environment2$pool_riffle),
as.numeric(environment2$meander))
colnames(X) <- c("avcondition", "avlength", "T", "O2", "Con", "COD", "NH4", "Nt", "netcen", "updist", "pool_riffle", "meander")
example_mcmc_control <- list(n.burnin = 1000, n.iteration = 10000, n.thin = 1)
testpath <- file.path(tempdir(), "jagsboralmodel.txt")
paramod <- boral(y, X = X,
family = "negative.binomial",
mcmc.control = example_mcmc_control,
model.name = testpath,
lv.control = list(num.lv = 2, type = "independent"),
save.model = TRUE)## module glm loaded## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 185
## Unobserved stochastic nodes: 338
## Total graph size: 2173
##
## Initializing modelplot(paramod)## NULL
coefsplot(covname = "avcondition", object = paramod) #Condition
coefsplot(covname = "avlength", object = paramod) #Length
coefsplot(covname = "T", object = paramod) #Temperature
coefsplot(covname = "O2", object = paramod) #Oxygen
coefsplot(covname = "Con", object = paramod) #Conductivity
coefsplot(covname = "COD", object = paramod) #COD
coefsplot(covname = "NH4", object = paramod) #NH4
coefsplot(covname = "Nt", object = paramod) #Nt
coefsplot(covname = "netcen", object = paramod) #netcen
coefsplot(covname = "updist", object = paramod) #updist
coefsplot(covname = "pool_riffle", object = paramod) #poolriffle
coefsplot(covname = "meander", object = paramod) #meander
envcors <- get.enviro.cor(paramod)
rescors <- get.residual.cor(paramod)
corrplot(envcors$sig.cor, type = "lower", diag = FALSE, title = "Correlations due to covariates", mar = c(3,0.5,2,1), tl.srt = 45)
corrplot(rescors$sig.cor, type = "lower", diag = FALSE, title = "Residual correlations", mar = c(3,0.5,2,1), tl.srt = 45)
7.1 BORAL analysis for median infection intensities of parasites
y <- round(cbind(medin$Gyr, medin$Tri, medin$Glu, medin$Con, medin$Ang))
paramod <- boral(y, X = X,
family = "negative.binomial",
mcmc.control = example_mcmc_control,
model.name = testpath,
lv.control = list(num.lv = 2, type = "independent"),
save.model = TRUE)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 185
## Unobserved stochastic nodes: 338
## Total graph size: 2173
##
## Initializing modelplot(paramod)## NULL
coefsplot(covname = "avcondition", object = paramod) #Condition
coefsplot(covname = "avlength", object = paramod) #Length
coefsplot(covname = "T", object = paramod) #Temperature
coefsplot(covname = "O2", object = paramod) #Oxygen
coefsplot(covname = "Con", object = paramod) #Conductivity
coefsplot(covname = "COD", object = paramod) #COD
coefsplot(covname = "NH4", object = paramod) #NH4
coefsplot(covname = "Nt", object = paramod) #Nt
coefsplot(covname = "netcen", object = paramod) #netcen
coefsplot(covname = "updist", object = paramod) #updist
coefsplot(covname = "pool_riffle", object = paramod) #poolriffle
coefsplot(covname = "meander", object = paramod) #meander
envcors <- get.enviro.cor(paramod)
rescors <- get.residual.cor(paramod)
corrplot(envcors$sig.cor, type = "lower", diag = FALSE, title = "Correlations due to covariates", mar = c(3,0.5,2,1), tl.srt = 45)
corrplot(rescors$sig.cor, type = "lower", diag = FALSE, title = "Residual correlations", mar = c(3,0.5,2,1), tl.srt = 45)
8. Multivariate analysis
8.1 Component communities
# Component communities: Bray-Curtis dissimilarities based on Hellinger transformed average abundance data
spe.hel_bray <- vegdist(decostand(avab[,-1], na.rm=T, method="hellinger"), method="bray", na.rm=T)
# Check whether Euclidean and Bray-Curtis distances are comparable
spe.hel_euc <- vegdist(decostand(avab[,-1], na.rm=T, method="hellinger"), method="euc", na.rm=T)
plot(spe.hel_bray, spe.hel_euc)
mantel(spe.hel_bray, spe.hel_euc)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = spe.hel_bray, ydis = spe.hel_euc)
##
## Mantel statistic r: 0.9648
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.110 0.144 0.186 0.206
## Permutation: free
## Number of permutations: 999model_adonis = adonis(spe.hel_bray ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2)## 'adonis' will be deprecated: use 'adonis2' insteadmodel_adonis$aov.tab## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## avlength 1 0.1863 0.18629 1.8398 0.04303 0.110
## avcondition 1 0.1305 0.13055 1.2893 0.03016 0.272
## T_av 1 0.3204 0.32039 3.1643 0.07401 0.020 *
## O2_sat_av 1 0.1368 0.13678 1.3509 0.03160 0.266
## Con_av 1 0.1326 0.13262 1.3098 0.03064 0.280
## COD_av 1 0.0657 0.06568 0.6486 0.01517 0.658
## NH4._av 1 0.2040 0.20399 2.0146 0.04712 0.110
## Nt_av 1 0.1451 0.14509 1.4329 0.03352 0.212
## pool_riffle 1 0.0853 0.08529 0.8424 0.01970 0.496
## meander 1 0.2029 0.20286 2.0035 0.04686 0.098 .
## netcen 1 0.2173 0.21728 2.1459 0.05019 0.071 .
## updist 1 0.0719 0.07193 0.7104 0.01662 0.577
## Residuals 24 2.4301 0.10125 0.56137
## Total 36 4.3288 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# environmental variables
env_select <- environment2[,c("T_av", "O2_sat_av", "Con_av", "COD_av", "NH4._av", "Nt_av", "pool_riffle", "meander", "netcen", "updist")]
env_select$pool_riffle <- as.numeric(env_select$pool_riffle)
env_select$meander <- as.numeric(env_select$meander)
pca <- prcomp(env_select, scale.=T)
summary(pca)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.7124 1.5545 1.1221 1.0140 0.88807 0.79463 0.56647
## Proportion of Variance 0.2933 0.2416 0.1259 0.1028 0.07887 0.06314 0.03209
## Cumulative Proportion 0.2933 0.5349 0.6608 0.7636 0.84248 0.90563 0.93771
## PC8 PC9 PC10
## Standard deviation 0.50483 0.46939 0.38429
## Proportion of Variance 0.02549 0.02203 0.01477
## Cumulative Proportion 0.96320 0.98523 1.00000plot(pca)
biplot(pca)
# Assess the effect of environmental variables on parasite component community dissimilarities using distance based RDA
spe.rda <- dbrda(spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander, data = env_select)
anova(spe.rda)## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 8 1.1664 1.2909 0.163
## Residual 28 3.1624RsquareAdj(spe.rda)$adj.r.squared## [1] 0.06072755mod0 <- dbrda(spe.hel_bray ~ 1, env_select[,-c(9:10)]) # Model with intercept only #edit_PH
mod1 <- dbrda(spe.hel_bray ~ ., env_select[,-c(9:10)]) # Model with all explanatory variables #edit_PH
step.res <- ordiR2step(mod0, mod1, direction = "both",perm.max = 200)## Step: R2.adj= 0
## Call: spe.hel_bray ~ 1
##
## R2.adjusted
## <All variables> 0.060727547
## + T_av 0.036398792
## + NH4._av 0.020208612
## + meander 0.018502880
## + O2_sat_av 0.004277611
## + Con_av 0.001872668
## + pool_riffle 0.001742860
## <none> 0.000000000
## + Nt_av -0.002060170
## + COD_av -0.017968936
##
## Df AIC F Pr(>F)
## + T_av 1 54.788 2.3599 0.048 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.03639879
## Call: spe.hel_bray ~ T_av
##
## R2.adjusted
## + NH4._av 0.07975512
## <All variables> 0.06072755
## + meander 0.05721602
## + Con_av 0.04792452
## + Nt_av 0.04664053
## + O2_sat_av 0.04174670
## + pool_riffle 0.03938572
## <none> 0.03639879
## + COD_av 0.02324854step.res$anova # Summary table## R2.adj Df AIC F Pr(>F)
## + T_av 0.036399 1 54.788 2.3599 0.048 *
## <All variables> 0.060728
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(spe.rda, scaling = 1) # it is for technical reasons not possible to plot both site and species scores
summary(spe.rda)##
## Call:
## dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
##
## Partitioning of squared Bray distance:
## Inertia Proportion
## Total 4.329 1.0000
## Constrained 1.166 0.2695
## Unconstrained 3.162 0.7305
##
## Eigenvalues, and their contribution to the squared Bray distance
##
## Importance of components:
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6
## Eigenvalue 0.5740 0.3372 0.17566 0.09415 0.042525 0.0020670
## Proportion Explained 0.1326 0.0779 0.04058 0.02175 0.009824 0.0004775
## Cumulative Proportion NA NA NA NA NA NA
## idbRDA1 idbRDA2 MDS1 MDS2 MDS3 MDS4 MDS5
## Eigenvalue -0.014797 -0.04441 1.3249 0.8027 0.4680 0.31323 0.29121
## Proportion Explained 0.003418 0.01026 0.3061 0.1854 0.1081 0.07236 0.06727
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS6 MDS7 MDS8 MDS9 MDS10 MDS11 MDS12
## Eigenvalue 0.13644 0.11438 0.08987 0.07065 0.06425 0.02766 0.01576
## Proportion Explained 0.03152 0.02642 0.02076 0.01632 0.01484 0.00639 0.00364
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS13 MDS14 iMDS1 iMDS2 iMDS3
## Eigenvalue 0.011872 2.725e-04 -0.0020314 -0.007356 -0.01078
## Proportion Explained 0.002743 6.295e-05 0.0004693 0.001699 0.00249
## Cumulative Proportion NA NA NA NA NA
## iMDS4 iMDS5 iMDS6 iMDS7 iMDS8 iMDS9
## Eigenvalue -0.014744 -0.01935 -0.022495 -0.029849 -0.03342 -0.041489
## Proportion Explained 0.003406 0.00447 0.005197 0.006895 0.00772 0.009584
## Cumulative Proportion NA NA NA NA NA NA
## iMDS10 iMDS11 iMDS12 iMDS13 iMDS14
## Eigenvalue -0.05984 -0.06439 -0.06995 -0.09116 -0.10203
## Proportion Explained 0.01382 0.01487 0.01616 0.02106 0.02357
## Cumulative Proportion NA NA NA NA NA
##
## Accumulated constrained eigenvalues
## Importance of components:
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6 idbRDA1
## Eigenvalue 0.5740 0.3372 0.1757 0.09415 0.04252 0.002067 -0.01480
## Proportion Explained 0.4921 0.2891 0.1506 0.08072 0.03646 0.001772 0.01269
## Cumulative Proportion NA NA NA NA NA NA NA
## idbRDA2
## Eigenvalue -0.04441
## Proportion Explained 0.03807
## Cumulative Proportion NA
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 3.533199
##
##
## Site scores (weighted sums of species scores)
##
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6
## 1 -0.07530 -0.62580 -0.63602 -2.643196 -0.34427 -23.8946
## 2 0.87242 0.77832 -0.32068 2.074749 -0.25800 41.3754
## 3 -0.81094 1.29200 -0.57735 -0.708420 -2.35800 -5.0035
## 4 -0.97450 0.91051 1.79849 -2.192492 -1.99248 -4.7134
## 5 -0.25889 0.74559 -1.30995 -1.418227 0.33185 -21.9720
## 6 0.41640 1.01888 -0.34551 -2.318324 1.61924 -24.8949
## 7 -0.98286 -0.67495 0.62017 1.174218 1.55155 -0.8974
## 8 0.93807 0.62635 -3.10878 -0.765073 2.41981 -3.0764
## 9 -0.62816 -0.77245 -1.30540 1.806170 1.30543 -41.2679
## 10 -1.03727 -0.50089 0.62572 0.369727 1.00378 4.7819
## 11 1.10743 -0.99448 -0.66515 -0.071830 0.83402 -14.8528
## 12 0.61624 0.76828 0.87110 0.357836 0.12125 38.2552
## 13 -0.96088 -0.72453 0.47977 1.569869 1.43399 -6.8692
## 14 1.41128 2.03453 -1.67256 0.776255 0.27648 29.9435
## 15 -1.00448 -0.59302 0.65010 0.833950 1.44643 3.4902
## 16 1.26382 0.81623 3.19230 1.380128 1.12523 -10.4959
## 17 -0.92726 -0.73637 0.33774 1.892276 1.07787 -8.2400
## 18 -1.04523 -1.01280 -0.59052 2.897871 0.35870 -14.2865
## 19 0.89987 -1.54255 -0.54435 1.738188 -2.61939 32.8867
## 20 -0.14608 -0.34079 0.78630 1.316880 1.70249 -9.1246
## 21 -1.30124 -1.46318 -1.27370 3.739798 0.91928 -16.7232
## 22 -0.07485 1.73636 0.30728 -1.760677 -0.46041 6.1504
## 23 -0.18008 -0.64105 0.55428 -1.310793 -0.82437 -3.6579
## 24 0.53054 0.08088 1.31623 -3.529128 -1.96164 41.6937
## 25 -0.50901 -0.82222 -1.26480 2.167458 -6.71178 -22.4535
## 26 0.54874 0.92050 1.12409 0.009258 0.07714 34.0773
## 27 1.25804 -0.69768 2.37255 -1.697654 -3.24936 -16.2349
## 28 -0.86554 0.32492 -0.68347 0.704689 -0.16329 13.0277
## 29 1.82276 0.98120 0.33891 -1.852448 4.33581 -8.6587
## 30 0.72974 -0.59907 -1.00943 -0.201449 1.77124 -4.5186
## 31 -1.02904 0.19515 1.05102 -1.468951 -1.36610 23.1359
## 32 -0.13877 -0.06050 -0.19569 -1.899083 -1.10403 -31.4856
## 33 1.01209 -0.59332 2.58400 -0.318337 -1.43566 39.0635
## 34 -0.39319 -1.79743 0.20497 -0.362972 -0.59693 -19.8004
## 35 -0.89506 1.46497 -0.05622 -1.876722 -1.86120 3.5844
## 36 0.87386 1.17335 -3.50567 1.196922 1.51036 -9.4970
## 37 -0.06267 -0.67495 -0.14975 0.389535 2.08495 11.1529
##
##
## Site constraints (linear combinations of constraining variables)
##
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6
## 1 -0.16605 -0.69736 -0.45052 0.39949 -1.140870 0.26766
## 2 -0.07683 -0.33692 0.61638 0.50606 -0.381728 1.66676
## 3 -0.10071 0.76299 0.51759 0.02910 0.003366 -0.26259
## 4 -1.55548 0.90013 0.66933 -0.51241 -0.415318 -0.13093
## 5 -0.41914 0.27994 -0.66104 -1.10909 0.384267 0.50060
## 6 0.63090 0.30135 -0.07724 0.37780 0.857007 -0.33310
## 7 -0.15252 -0.25647 -0.20458 -0.38652 0.792931 0.48033
## 8 -0.14874 0.12938 -0.39866 -0.78891 0.083814 0.27677
## 9 -0.14359 0.11628 0.95439 -0.04629 0.609992 -0.64671
## 10 -0.74154 -0.65088 0.87072 -1.17151 -0.875610 -1.73263
## 11 0.77571 -0.48314 0.02321 -0.76003 0.259283 -0.01494
## 12 -0.44387 -0.28442 0.54938 -0.01813 1.062736 -0.71676
## 13 -0.25608 -0.61494 -0.47847 0.74269 0.180828 -0.14346
## 14 -0.23636 0.67472 -0.87185 0.72759 0.046993 -0.61306
## 15 -0.75508 0.19781 0.41035 0.89449 -0.744271 -0.25040
## 16 0.97946 1.50039 0.74911 0.12808 -0.657115 -0.30363
## 17 -1.19679 -0.18091 -0.90174 0.16093 -0.229189 0.70454
## 18 0.23197 0.39183 -0.42932 -0.82808 -0.242760 0.64991
## 19 0.80868 -0.55920 -0.22435 0.32426 -0.445610 -0.51746
## 20 -0.19206 0.12003 0.24728 0.99541 0.747110 0.40576
## 21 -0.38477 0.09371 -0.03632 1.07969 1.383306 -0.25819
## 22 0.02352 0.34014 0.07462 0.03662 -0.352806 0.35626
## 23 0.18716 -0.91860 -0.02126 -0.69548 0.117329 0.60753
## 24 0.03851 0.57991 -0.70662 -1.10691 0.533700 -0.10095
## 25 0.61401 -0.01017 -0.51875 0.25795 -1.076732 -0.24901
## 26 0.42275 0.07490 0.99168 0.37950 -0.249480 0.19511
## 27 0.54478 -0.12356 0.31406 0.41569 0.237094 0.09772
## 28 0.02044 0.43106 -0.71781 0.33274 -0.428303 -0.26278
## 29 1.23235 0.15293 0.10463 -0.07285 0.107338 0.07885
## 30 0.02214 -1.33353 -0.37271 0.38762 -0.062960 -0.48751
## 31 -0.71637 0.11204 -0.85950 0.24734 -0.066602 -0.17455
## 32 0.01756 -0.71899 -0.27770 0.01727 0.106473 -0.73172
## 33 0.78296 -0.75806 0.27650 -0.27767 -0.031297 -0.01895
## 34 -0.33590 -0.63554 1.07989 -0.07860 -0.526645 0.97933
## 35 -0.15761 0.49842 0.66426 -0.05091 0.331527 0.71522
## 36 0.53360 0.88375 -0.84613 0.03571 -0.630019 -0.08456
## 37 0.31300 0.02099 -0.05877 -0.57261 0.712218 0.05155
##
##
## Biplot scores for constraining variables
##
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5 dbRDA6
## T_av 0.63982 0.25527 0.1662 -0.34674 -0.2289 -0.48193
## O2_sat_av 0.02371 -0.02938 -0.8854 -0.18557 0.2035 -0.13356
## Con_av -0.25315 -0.27391 -0.1407 -0.76998 -0.4264 0.07657
## COD_av -0.13235 0.15034 0.4807 0.08407 0.1430 -0.57341
## NH4._av -0.40051 0.14159 0.7569 -0.21438 -0.2051 -0.37005
## Nt_av -0.31714 -0.23687 0.2952 -0.33830 -0.5347 -0.40917
## pool_riffle -0.34609 -0.15029 -0.4384 0.27287 -0.5547 -0.42271
## meander -0.33563 0.52703 -0.3800 -0.21165 -0.5422 -0.13159anova(spe.rda)## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 8 1.1664 1.2909 0.167
## Residual 28 3.1624anova(spe.rda, by="term")## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
## Df SumOfSqs F Pr(>F)
## T_av 1 0.2734 2.4210 0.048 *
## O2_sat_av 1 0.1377 1.2195 0.302
## Con_av 1 0.1613 1.4283 0.237
## COD_av 1 0.0676 0.5990 0.695
## NH4._av 1 0.1709 1.5131 0.212
## Nt_av 1 0.0501 0.4439 0.787
## pool_riffle 1 0.0657 0.5818 0.701
## meander 1 0.2396 2.1210 0.070 .
## Residual 28 3.1624
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova.cca(spe.rda, step=1000);## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 8 1.1664 1.2909 0.164
## Residual 28 3.1624anova.cca(spe.rda, step=1000, by="term");## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
## Df SumOfSqs F Pr(>F)
## T_av 1 0.2734 2.4210 0.057 .
## O2_sat_av 1 0.1377 1.2195 0.314
## Con_av 1 0.1613 1.4283 0.245
## COD_av 1 0.0676 0.5990 0.705
## NH4._av 1 0.1709 1.5131 0.218
## Nt_av 1 0.0501 0.4439 0.780
## pool_riffle 1 0.0657 0.5818 0.699
## meander 1 0.2396 2.1210 0.086 .
## Residual 28 3.1624
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1RsquareAdj(spe.rda)$adj.r.squared;## [1] 0.06072755RsquareAdj(spe.rda)$r.squared## [1] 0.26945488.1.2 Effect of space on component community structure
# Same for spatial predictors
spe.rda <- dbrda(spe.hel_bray ~ netcen + updist, data = env_select)
anova(spe.rda)## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ netcen + updist, data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 2 0.5154 2.2975 0.036 *
## Residual 34 3.8135
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1RsquareAdj(spe.rda)$adj.r.squared## [1] 0.06723433mod0 <- dbrda(spe.hel_bray ~ 1, env_select[,c(9:10)]) # Model with intercept only #edit_PH
mod1 <- dbrda(spe.hel_bray ~ ., env_select[,c(9:10)]) # Model with all explanatory variables #edit_PH
step.res <- ordiR2step(mod0, mod1, direction = "both",perm.max = 200)## Step: R2.adj= 0
## Call: spe.hel_bray ~ 1
##
## R2.adjusted
## <All variables> 0.06723433
## + updist 0.04867127
## + netcen 0.04317133
## <none> 0.00000000
##
## Df AIC F Pr(>F)
## + updist 1 54.314 2.8418 0.024 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.04867127
## Call: spe.hel_bray ~ updist
##
## R2.adjusted
## <All variables> 0.06723433
## + netcen 0.06723433
## <none> 0.04867127
##
## Df AIC F Pr(>F)
## + netcen 1 54.512 1.6965 0.152step.res$anova # Summary table## R2.adj Df AIC F Pr(>F)
## + updist 0.048671 1 54.314 2.8418 0.024 *
## <All variables> 0.067234
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1RsquareAdj(spe.rda)$adj.r.squared## [1] 0.06723433anova.cca(spe.rda, step=1000);## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ netcen + updist, data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 2 0.5154 2.2975 0.023 *
## Residual 34 3.8135
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova.cca(spe.rda, step=1000, by="term");## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ netcen + updist, data = env_select)
## Df SumOfSqs F Pr(>F)
## netcen 1 0.3019 2.6920 0.033 *
## updist 1 0.2134 1.9029 0.113
## Residual 34 3.8135
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1RsquareAdj(spe.rda)$adj.r.squared;## [1] 0.06723433RsquareAdj(spe.rda)$r.squared## [1] 0.11905468.1.3 Variation partitioning
#Variation partitioning
spe.varpart1 <- varpart(spe.hel_bray, env_select[,1:8], env_select[,9:10])
plot(spe.varpart1,digits=2)
spe.varpart1##
## Partition of squared Bray distance in dbRDA
##
## Call: varpart(Y = spe.hel_bray, X = env_select[, 1:8], env_select[,
## 9:10])
##
## Explanatory tables:
## X1: env_select[, 1:8]
## X2: env_select[, 9:10]
##
## No. of explanatory tables: 2
## Total variation (SS): 4.3288
## No. of observations: 37
##
## Partition table:
## Df R.squared Adj.R.squared Testable
## [a+c] = X1 8 0.26945 0.06073 TRUE
## [b+c] = X2 2 0.11905 0.06723 TRUE
## [a+b+c] = X1+X2 10 0.36716 0.12376 TRUE
## Individual fractions
## [a] = X1|X2 8 0.05653 TRUE
## [b] = X2|X1 2 0.06303 TRUE
## [c] 0 0.00420 FALSE
## [d] = Residuals 0.87624 FALSE
## ---
## Use function 'dbrda' to test significance of fractions of interestanova.cca(dbrda(spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + Condition(netcen + updist),
data=env_select), step=1000)## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander + Condition(netcen + updist), data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 8 1.0740 1.2742 0.195
## Residual 26 2.7395anova.cca(dbrda(spe.hel_bray ~ netcen + updist+
Condition(T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander), data=env_select), step=1000)## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = spe.hel_bray ~ netcen + updist + Condition(T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander), data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 2 0.42295 2.0071 0.047 *
## Residual 26 2.73945
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 18.2 Infra-communities
# Infracommunities: Bray-Curtis dissimilarities are calculated at the individual host level Hellinger-transformed parasite data and then averaged within site
# A dummy parasite species is added to avoid problems with non-infected fishes
data_infra <- na.omit(data[,c(1,22:24,26:32)])
data_infra_disp <- dispweight(data_infra[,-1])
braycurtis <- vegdist(decostand(cbind(data_infra_disp,rep(1,nrow(data_infra))), na.rm=T, method="hellinger"), method="bray", na.rm=T)
meandist_bray <- meandist(braycurtis, data_infra[,1])
# Check whether Euclidean and Bray-Curtis distances are comparable
braycurtis <- vegdist(decostand(cbind(data_infra_disp,rep(1,nrow(data_infra))), na.rm=T, method="hellinger"), method="bray", na.rm=T)
meandist_bray <- meandist(braycurtis, data_infra[,1])
euc <- vegdist(decostand(cbind(data_infra_disp,rep(1,nrow(data_infra))), na.rm=T, method="hellinger"), method="euc", na.rm=T)
meandist_euc <- meandist(euc, data_infra[,1])
plot(meandist_bray[1:37,1:37], meandist_euc[1:37,1:37])
mantel(meandist_bray[1:37,1:37], meandist_euc[1:37,1:37])##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = meandist_bray[1:37, 1:37], ydis = meandist_euc[1:37, 1:37])
##
## Mantel statistic r: 0.9906
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.184 0.241 0.292 0.391
## Permutation: free
## Number of permutations: 999model_adonis = adonis(meandist_bray ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + netcen +
updist, data=environment2)## 'adonis' will be deprecated: use 'adonis2' insteadmodel_adonis$aov.tab## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## avlength 1 0.05123 0.051226 6.1160 0.11570 0.006 **
## avcondition 1 0.00393 0.003927 0.4688 0.00887 0.622
## T_av 1 0.00874 0.008741 1.0437 0.01974 0.369
## O2_sat_av 1 0.02116 0.021158 2.5262 0.04779 0.105
## Con_av 1 0.05861 0.058606 6.9971 0.13237 0.006 **
## COD_av 1 0.03543 0.035428 4.2298 0.08002 0.027 *
## NH4._av 1 0.00275 0.002754 0.3288 0.00622 0.690
## Nt_av 1 0.01049 0.010488 1.2521 0.02369 0.293
## pool_riffle 1 0.00491 0.004914 0.5867 0.01110 0.563
## meander 1 0.00959 0.009588 1.1447 0.02166 0.347
## netcen 1 0.02319 0.023187 2.7684 0.05237 0.070 .
## updist 1 0.01171 0.011707 1.3977 0.02644 0.283
## Residuals 24 0.20102 0.008376 0.45403
## Total 36 0.44274 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# environmental variables
env_select <- environment2[,c("T_av", "O2_sat_av", "Con_av", "COD_av", "NH4._av", "Nt_av", "pool_riffle", "meander", "netcen", "updist")]
env_select$pool_riffle <- as.numeric(env_select$pool_riffle)
env_select$meander <- as.numeric(env_select$meander)
pca <- prcomp(env_select, scale.=T)
summary(pca)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.7124 1.5545 1.1221 1.0140 0.88807 0.79463 0.56647
## Proportion of Variance 0.2933 0.2416 0.1259 0.1028 0.07887 0.06314 0.03209
## Cumulative Proportion 0.2933 0.5349 0.6608 0.7636 0.84248 0.90563 0.93771
## PC8 PC9 PC10
## Standard deviation 0.50483 0.46939 0.38429
## Proportion of Variance 0.02549 0.02203 0.01477
## Cumulative Proportion 0.96320 0.98523 1.00000plot(pca)
biplot(pca)
8.2.1 Effect of environment on infracommunity structure
# Assess the effect of environmental variables on parasite infracommunity dissimilarities using distance based RDA
spe.rda <- dbrda(meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander, data = env_select)
anova(spe.rda)## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 8 0.45036 1.3833 0.008 **
## Residual 28 1.13946
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1RsquareAdj(spe.rda)$adj.r.squared## [1] 0.07850103mod0 <- dbrda(meandist_bray ~ 1, env_select) # Model with intercept only #edit_PH
mod1 <- dbrda(meandist_bray ~ ., env_select) # Model with all explanatory variables #edit_PH
step.res <- ordiR2step(mod0, mod1, direction = "both",perm.max = 200)## Step: R2.adj= 0
## Call: meandist_bray ~ 1
##
## R2.adjusted
## <All variables> 0.0898822750
## + Con_av 0.0492715387
## + NH4._av 0.0376899022
## + Nt_av 0.0353578268
## + COD_av 0.0097867139
## + updist 0.0092890919
## + pool_riffle 0.0070398050
## + netcen 0.0034499960
## + O2_sat_av 0.0031240320
## <none> 0.0000000000
## + T_av -0.0005031246
## + meander -0.0037253892
##
## Df AIC F Pr(>F)
## + Con_av 1 17.228 2.8657 0.004 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.04927154
## Call: meandist_bray ~ Con_av
##
## R2.adjusted
## <All variables> 0.08988227
## + COD_av 0.08426096
## + NH4._av 0.07058808
## + updist 0.06800179
## + O2_sat_av 0.06266812
## + netcen 0.05455352
## + meander 0.05428744
## + Nt_av 0.05273614
## <none> 0.04927154
## + pool_riffle 0.04923534
## + T_av 0.04694754
##
## Df AIC F Pr(>F)
## + COD_av 1 16.768 2.3373 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.08426096
## Call: meandist_bray ~ Con_av + COD_av
##
## R2.adjusted
## + updist 0.09961037
## + meander 0.09248673
## <All variables> 0.08988227
## + netcen 0.08626581
## + pool_riffle 0.08559576
## <none> 0.08426096
## + T_av 0.08115910
## + O2_sat_av 0.08092067
## + Nt_av 0.07979657
## + NH4._av 0.07591794step.res$anova # Summary table## R2.adj Df AIC F Pr(>F)
## + Con_av 0.049272 1 17.228 2.8657 0.004 **
## + COD_av 0.084261 1 16.768 2.3373 0.002 **
## <All variables> 0.089882
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1spe.rda <- dbrda(meandist_bray ~ Con_av + COD_av, env_select)
plot(spe.rda, scaling = 1) # it is for technical reasons not possible to plot both site and species scores
summary(spe.rda)##
## Call:
## dbrda(formula = meandist_bray ~ Con_av + COD_av, data = env_select)
##
## Partitioning of squared Unknown distance:
## Inertia Proportion
## Total 1.5898 1.0000
## Constrained 0.2148 0.1351
## Unconstrained 1.3750 0.8649
##
## Eigenvalues, and their contribution to the squared Unknown distance
##
## Importance of components:
## dbRDA1 dbRDA2 MDS1 MDS2 MDS3 MDS4 MDS5
## Eigenvalue 0.1863 0.02858 0.1972 0.15873 0.1207 0.08661 0.07076
## Proportion Explained 0.1172 0.01798 0.1240 0.09984 0.0759 0.05447 0.04451
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS6 MDS7 MDS8 MDS9 MDS10 MDS11 MDS12
## Eigenvalue 0.06231 0.04802 0.04613 0.04387 0.04231 0.03826 0.03703
## Proportion Explained 0.03919 0.03021 0.02901 0.02760 0.02661 0.02406 0.02329
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS13 MDS14 MDS15 MDS16 MDS17 MDS18 MDS19
## Eigenvalue 0.03557 0.03355 0.03200 0.02808 0.02797 0.02575 0.02487
## Proportion Explained 0.02237 0.02110 0.02013 0.01766 0.01759 0.01620 0.01564
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS20 MDS21 MDS22 MDS23 MDS24 MDS25 MDS26
## Eigenvalue 0.02377 0.02238 0.02182 0.02104 0.01990 0.01880 0.01782
## Proportion Explained 0.01495 0.01408 0.01372 0.01323 0.01252 0.01183 0.01121
## Cumulative Proportion NA NA NA NA NA NA NA
## MDS27 MDS28 MDS29 MDS30 MDS31 MDS32
## Eigenvalue 0.01719 0.01622 0.011778 0.010499 0.007903 0.006958
## Proportion Explained 0.01081 0.01020 0.007408 0.006604 0.004971 0.004376
## Cumulative Proportion NA NA NA NA NA NA
## MDS33 iMDS1
## Eigenvalue 0.005670 -0.006426
## Proportion Explained 0.003566 0.004042
## Cumulative Proportion NA NA
##
## Accumulated constrained eigenvalues
## Importance of components:
## dbRDA1 dbRDA2
## Eigenvalue 0.1863 0.02858
## Proportion Explained 0.8670 0.13304
## Cumulative Proportion 0.8670 1.00000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 2.750505
##
##
## Site scores (weighted sums of species scores)
##
## dbRDA1 dbRDA2 MDS1 MDS2 MDS3 MDS4
## SITE 1 -0.215744 0.44381 0.01302 0.31819 -0.28116 -0.178171
## SITE 11 -0.507248 0.82144 -0.09719 -0.28198 -0.04679 0.342258
## SITE 12 -0.746029 0.15775 0.31672 0.25334 -0.08234 -0.151368
## SITE 13 -0.005847 -0.21986 -0.02027 0.91847 -0.12233 -0.332031
## SITE 14 -0.101987 0.28502 0.01314 0.55873 -0.09066 -0.223950
## SITE 15 -0.290401 -0.51714 -0.23614 -0.10408 -0.36670 -0.263022
## SITE 16 0.016467 0.22318 0.14693 -0.09879 -0.21010 0.156923
## SITE 17 -0.053594 0.78979 0.68552 0.22623 0.55763 0.068555
## SITE 18 0.741297 -0.34140 0.34156 -0.78781 -0.33039 0.372025
## SITE 19 2.204147 0.31903 0.07576 -0.77176 -0.32619 0.237586
## SITE 2 -0.334822 0.53901 0.20181 0.23760 0.25255 0.460075
## SITE 20 0.132501 -1.43824 -0.36307 -0.15683 0.91526 0.652586
## SITE 21 -0.518164 0.15765 0.51672 0.22690 -0.05180 0.071838
## SITE 22 -0.681351 -0.71637 0.44536 0.61813 0.24285 -0.440713
## SITE 23 0.085803 -0.07384 -0.18742 -0.58764 -0.50280 0.146049
## SITE 24 -0.506368 -1.41199 -0.02580 0.53865 0.25884 0.040947
## SITE 26 0.017028 0.61108 0.37526 -0.36715 -0.10492 0.425754
## SITE 28 0.121498 0.47609 0.59944 -0.22080 -0.07107 0.465237
## SITE 29 -0.014961 0.35365 0.04914 -0.10076 -0.08378 0.641604
## SITE 3 -0.711132 -0.09645 0.08469 -0.08718 -0.21517 0.038217
## SITE 30 -0.812113 -0.12175 0.55582 -0.16918 -0.25866 0.087617
## SITE 31 0.043616 -0.61571 -0.31459 0.47663 -0.16031 -0.651535
## SITE 32 0.562390 0.35603 -0.38737 0.06837 -0.20621 -0.072343
## SITE 33 -0.372173 0.03501 0.10319 0.85337 -0.01445 -0.328696
## SITE 34 -0.466651 0.48378 0.50738 0.02349 -0.08407 0.306511
## SITE 35 0.163917 -1.66724 -0.04520 0.22652 0.57157 0.198776
## SITE 36 -0.296043 -0.76361 -0.59812 0.40106 -0.04152 -0.249681
## SITE 38 0.302702 0.34509 -0.16274 -0.49724 -0.25040 -0.005416
## SITE 39 -0.978171 0.31142 0.36167 0.40699 0.02355 0.020399
## SITE 4 -0.141634 0.44343 -0.02774 -0.11080 -0.07834 0.180628
## SITE 40 0.496377 -0.01681 -0.38350 -0.04028 -0.50240 -0.167352
## SITE 41 1.117637 0.21808 -0.82799 -0.60061 -0.49878 -0.637014
## SITE 42 0.437517 -0.97421 -1.46550 0.08641 1.33505 0.880493
## SITE 5 0.522987 0.17736 0.27513 0.18757 -0.04206 0.152978
## SITE 6 0.423826 -0.07352 -0.91744 -0.10155 -0.42928 -0.789691
## SITE 7 0.013170 1.17661 0.57227 -1.15454 1.52741 -1.616432
## SITE 9 0.351552 0.32386 -0.18044 -0.38769 -0.23204 0.160360
##
##
## Site constraints (linear combinations of constraining variables)
##
## dbRDA1 dbRDA2 MDS1 MDS2 MDS3 MDS4
## SITE 1 -0.079483 0.44275 0.01302 0.31819 -0.28116 -0.178171
## SITE 11 -0.692558 0.69501 -0.09719 -0.28198 -0.04679 0.342258
## SITE 12 -0.387516 -0.29270 0.31672 0.25334 -0.08234 -0.151368
## SITE 13 0.692033 0.28074 -0.02027 0.91847 -0.12233 -0.332031
## SITE 14 0.273072 0.58517 0.01314 0.55873 -0.09066 -0.223950
## SITE 15 -0.503973 -0.49304 -0.23614 -0.10408 -0.36670 -0.263022
## SITE 16 -0.041630 0.31054 0.14693 -0.09879 -0.21010 0.156923
## SITE 17 0.286975 0.22462 0.68552 0.22623 0.55763 0.068555
## SITE 18 0.355899 -0.83587 0.34156 -0.78781 -0.33039 0.372025
## SITE 19 1.702755 -0.40556 0.07576 -0.77176 -0.32619 0.237586
## SITE 2 -0.068404 0.54101 0.20181 0.23760 0.25255 0.460075
## SITE 20 0.005618 -0.62143 -0.36307 -0.15683 0.91526 0.652586
## SITE 21 -0.085484 -0.08616 0.51672 0.22690 -0.05180 0.071838
## SITE 22 -0.015359 -0.63175 0.44536 0.61813 0.24285 -0.440713
## SITE 23 -0.475351 -0.06131 -0.18742 -0.58764 -0.50280 0.146049
## SITE 24 -0.138656 -0.63048 -0.02580 0.53865 0.25884 0.040947
## SITE 26 -0.086478 0.70999 0.37526 -0.36715 -0.10492 0.425754
## SITE 28 0.227529 0.47771 0.59944 -0.22080 -0.07107 0.465237
## SITE 29 -0.053389 0.16940 0.04914 -0.10076 -0.08378 0.641604
## SITE 3 -0.805011 -0.44279 0.08469 -0.08718 -0.21517 0.038217
## SITE 30 -0.711357 -0.97402 0.55582 -0.16918 -0.25866 0.087617
## SITE 31 0.197349 -0.41392 -0.31459 0.47663 -0.16031 -0.651535
## SITE 32 0.342870 0.50733 -0.38737 0.06837 -0.20621 -0.072343
## SITE 33 0.260419 0.25833 0.10319 0.85337 -0.01445 -0.328696
## SITE 34 -0.213627 0.34583 0.50738 0.02349 -0.08407 0.306511
## SITE 35 0.299656 -0.71657 -0.04520 0.22652 0.57157 0.198776
## SITE 36 -0.265626 -0.39128 -0.59812 0.40106 -0.04152 -0.249681
## SITE 38 -0.177490 0.11612 -0.16274 -0.49724 -0.25040 -0.005416
## SITE 39 -0.502932 0.21901 0.36167 0.40699 0.02355 0.020399
## SITE 4 -0.213933 0.31087 -0.02774 -0.11080 -0.07834 0.180628
## SITE 40 0.157625 -0.02847 -0.38350 -0.04028 -0.50240 -0.167352
## SITE 41 0.308647 0.03805 -0.82799 -0.60061 -0.49878 -0.637014
## SITE 42 -0.043511 0.13560 -1.46550 0.08641 1.33505 0.880493
## SITE 5 0.774214 -0.17770 0.27513 0.18757 -0.04206 0.152978
## SITE 6 -0.221478 0.23059 -0.91744 -0.10155 -0.42928 -0.789691
## SITE 7 -0.080948 0.35692 0.57227 -1.15454 1.52741 -1.616432
## SITE 9 -0.020466 0.24745 -0.18044 -0.38769 -0.23204 0.160360
##
##
## Biplot scores for constraining variables
##
## dbRDA1 dbRDA2 MDS1 MDS2 MDS3 MDS4
## Con_av 0.7628 0.6467 0 0 0 0
## COD_av 0.4413 -0.8974 0 0 0 0anova(spe.rda)## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = meandist_bray ~ Con_av + COD_av, data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 2 0.21484 2.6563 0.001 ***
## Residual 34 1.37498
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(spe.rda, by="term")## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = meandist_bray ~ Con_av + COD_av, data = env_select)
## Df SumOfSqs F Pr(>F)
## Con_av 1 0.12032 2.9752 0.001 ***
## COD_av 1 0.09452 2.3373 0.008 **
## Residual 34 1.37498
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 18.2.2 Effect of space on infracommunity structure
spe.rda <- dbrda(meandist_bray ~ netcen + updist, data = env_select)
anova(spe.rda)## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = meandist_bray ~ netcen + updist, data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 2 0.10742 1.2319 0.158
## Residual 34 1.48240RsquareAdj(spe.rda)$adj.r.squared## [1] 0.01271734anova.cca(spe.rda, step=1000);## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = meandist_bray ~ netcen + updist, data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 2 0.10742 1.2319 0.15
## Residual 34 1.48240anova.cca(spe.rda, step=1000, by="term");## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = meandist_bray ~ netcen + updist, data = env_select)
## Df SumOfSqs F Pr(>F)
## netcen 1 0.04949 1.1352 0.264
## updist 1 0.05792 1.3285 0.135
## Residual 34 1.48240RsquareAdj(spe.rda)$adj.r.squared;## [1] 0.01271734RsquareAdj(spe.rda)$r.squared## [1] 0.067566378.2.3 Variation partitioning
#Variation partitioning
spe.varpart1 <- varpart(meandist_bray, env_select[,1:8], env_select[,9:10])
plot(spe.varpart1,digits=2)
spe.varpart1##
## Partition of squared Unknown user-supplied distance in dbRDA
##
## Call: varpart(Y = meandist_bray, X = env_select[, 1:8], env_select[,
## 9:10])
##
## Explanatory tables:
## X1: env_select[, 1:8]
## X2: env_select[, 9:10]
##
## No. of explanatory tables: 2
## Total variation (SS): 1.5898
## No. of observations: 37
##
## Partition table:
## Df R.squared Adj.R.squared Testable
## [a+c] = X1 8 0.28328 0.07850 TRUE
## [b+c] = X2 2 0.06757 0.01272 TRUE
## [a+b+c] = X1+X2 10 0.34269 0.08988 TRUE
## Individual fractions
## [a] = X1|X2 8 0.07716 TRUE
## [b] = X2|X1 2 0.01138 TRUE
## [c] 0 0.00134 FALSE
## [d] = Residuals 0.91012 FALSE
## ---
## Use function 'dbrda' to test significance of fractions of interestanova.cca(dbrda(meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander + Condition(netcen + updist),
data=env_select), step=1000)## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander + Condition(netcen + updist), data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 8 0.4374 1.3603 0.014 *
## Residual 26 1.0450
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova.cca(dbrda(meandist_bray ~ netcen + updist+
Condition(T_av + O2_sat_av + Con_av + COD_av
+ NH4._av + Nt_av + pool_riffle + meander), data=env_select), step=1000)## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = meandist_bray ~ netcen + updist + Condition(T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander), data = env_select)
## Df SumOfSqs F Pr(>F)
## Model 2 0.09446 1.1751 0.175
## Residual 26 1.04500